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ABSTRACT 


The Simple Biosphere Model (SiB) of Sellers et al. (1986) was designed to simulate 
the interactions between the Earth’s land surface and the atmosphere by treating the 
vegetation explicitly and realistically, thereby incorporating biophysical controls on the 
exchanges of radiation, momentum, sensible and latent heat between the two systems. 
This paper describes the steps taken to implement SiB in a modified version of the 
National Meteorological Center’s spectral GCM. The coupled model (SiB-GCM) was used 
to produce summer and winter simulations. The same GCM was used with a conventional 
hydrological model (Ctl-GCM) to produce comparable ‘control’ summer and winter 
simulations. 

It was found that SiB-GCM produced a more realistic partitioning of energy at the 
land surface than Ctl-GCM. Generally, SiB-GCM produced more sensible heat flux and 
less latent heat flux over vegetated land than did Ctl-GCM and this resulted in the 
development of a much deeper daytime planetary boundary layer and reduced precipitation 
rates over the continents in SiB-GCM. In the summer simulation, the 200 mb jet stream 
and the wind speed at 850 mb were slightly weakened in the SiB-GCM relative to the 
Ctl-GCM results and equivalent analyses made from observations. 
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1. INTRODUCTION 

This report and the paper of Sato et al. (1989) describe a study in which a simple 
but biophysically based model of the terrestrial vegetation was implemented in a General 
Circulation Model (GCM) of the atmosphere. The motivation for the study was to assess 
the impact of implementing such a scheme on the simulated climate with particular 
reference to (i) observations where available and, (ii) a parallel simulation produced with 
an abiotic land surface model. The paper of Sato et al. (1989) reported on the results of 
the study: this report is complementary to the the paper in providing a detailed review of 
the technical aspects of the implementation and an expanded discussion of the reuslts. 

Over the last decade, the scientific community has recognized that the interaction 
between the land surface and the atmosphere may play a significant role in the climate 
system; specifically, the fluxes of radiation, momentum and sensible and latent heat 
between the two systems may influence the atmospheric motion, temperature, humidity 
and precipitation fields which in turn feed back onto the spatial distribution of the sinks 
and sources of these fluxes on the land surface. 

Until recently, the surface properties which regulate the exchange of radiation, 
momentum and heat between the land and the atmosphere have been regarded as separate 
items which may be independently prescribed as boundary conditions within a GCM. 
Although this approach is not very realistic, it lends itself to the investigation of land 
surface— atmosphere interactions by way of sensitivity studies conducted with GCM’s. The 
three land surface properties which ultimately govern these interactions are the albedo 
(radiative transfer), surface roughness (momentum transfer) and the surface hydrological 
parameterization/prescription (sensible and latent heat transfer), all of which have been 
the subject of individual sensitivity studies. 

First, the land surface albedo is known to vary from about 12% (tropical rain 
forest) to around 35% (Sahara desert) over the snowfree land, see Shuttleworth et al. 
(1984b) and Matthews (1984). Charney et al. (1977) executed a series of GCM sensitivity 
experiments centered on albedo changes in the Sahel which indicated that an increase in 
the albedo of the region would lead to a decrease in the surface evaporation rate and a 
reduction in the precipitation rate in the same area. This study has been followed by 
similar ones on the Sahel/Sahara region (Sud and Fennessy (1982) and Chervin (1979)) 
which produced the same result qualitatively if not quantitatively. Carson and Sangster 
(1981) performed a more radical experiment by comparing two simulation runs performed 
with a GCM wherein the global land surface albedo, a , was prescribed as 0.1 or 0.3. The 
higher albedo case was associated with a lowered precipitation rate over the land, (4.6 
mm/day for a = 0.1 and 3.4 mm/day for a = 0.3). 

The effect of changing surface roughness length (z 0 ) fields was investigated in a 
study by Sud and Smith (1985): The roughness length of the world’s deserts was reduced 
from 45 cm to 0.2 mm which brought about changes in the fields of horizontal water vapor 
convergence and convective precipitation. 

Not surprisingly, many sensitivity studies have been devoted to investigating the 
role of land surface hydrology in the exchange of energy between the surface and 
atmosphere. Walker and Rowntree (1977), Carson and Sangster (1981), Shukla and Mintz 
(1982), Warrilow (1986) and others reviewed in Mintz (1984) have performed experiments 
ranging from global prescriptions of totally wet versus totally dry land surfaces to 
alterations in regional soil moisture capacities. In almost all cases, reduced land surface 
evaporation rates led to reduced precipitation rates in the continental interiors. 



All of the above experiments served to demonstrate that the specification of albedo, 
roughness and perhaps most importantly, the ‘surface wetness’ could have important 
impacts on atmospheric fields. It was pointed out by Dickinson (1984) that these 
experiments, though illustrative, were unrealistic, particularly witn regard to the 
description of the land surface evapotranspiration process as conventionally specified in 
GCM’s. Dickinson (1984) designed a land surface model which explicitly refers to the 
effects of vegetation which has since been used for a simulation study of the climatic 
impact of large scale deforestation of the Amazon Basin; see Dickinson and Henderson— 
Sellers (1988). Sellers et al. (1986) constructed the Simple Biosphere model (SiB) which 
incorporates much of the philosophy of the Dickinson (1984) model while differing in detail. 
Both of these models have undergone sensitivity studies which indicate that as well as 
being realistic, the recognition of the role of vegetation in GCM’s should lead to more 
reasonable descriptions of the land surface— atmosphere fluxes. This paper reports on one 
such study, in which the conventional land surface hydrology model used in a research 
version of the National Meteorological Center (NMC) GCM was replaced by SiB. Changes 
in the simulated surface and atmospheric fields as a result of this implementation are 
discussed with reference to relevant observations where available. 

2. DESCRIPTION OF THE MODEL AND EXPERIMENTAL DESIGN 

The National Meteorological Center (NMC) spectral general circulation model 
(GCM) was originally designed and applied as a medium range forecast (MRF) model, see 
Sela (1980), but has recently been modified by us for use as a meteorological research tool. 
These modifications are reviewed below and include the implementation of a more efficient 
radiation scheme by Harshvardhan et al. (1987); the incorporation of a simple curve— fitted 
function in place of the full Monin— Obukhov scheme used to calculate surface transfer 
coefficients and resistances; inclusion of the Mellor and Yamada (1982) vertical diffusion 
scheme; and improvement of the time integration scheme used for the lower atmosphere 
and surface physics. 

The model formulation of dynamics is described in Sela (1980) and Kinter et al. 
(1988) describe the physics, initialization procedures and boundary conditions. The model 
has a rhomboidal truncation at total wave number 40, Laplacian-square type horizontal 
diffusion to avoid spectral blocking and is discretized into 18 vertical levels of which the 
PBL may occupy five or six. It has a shallow convection scheme (Tiedtke et al. (1984) as 
well as a large-scale precipitation and cumulus convection scheme after Kuo (1965) wnich 
serves as a moist process parameterization. The model has three soil layers with prognostic 
temperatures at 2.5 cm, 10 cm, and 50 cm depth; the temperature at 500 cm is kept 
constant. The surface hydrology is represented by a bucket model (Manabe, 1969) wherein 
the wetness is defined as the soil moisture divided by three quarters of maximum field 
capacity which is taken to be 150 mm on all land surfaces. These model ‘physics’ were 
taken from the 1986 version of NMC medium range forecast model. 

Some parameterizations were changed: 

1. A more efficient radiation scheme was incorporated which allowed simulation of the 
diurnal cycle; see Harshvardhan et al. (1987). The long wave radiation scheme is 
based on that of Harshvardhan and Corsetti (1984) and the short wave (solar) 
radiation scheme is based on that of Lads and Hansen (1974) as modified by Davies 
(1982). The solar radiation scheme was further modified to allow for partial 
cloudiness. The new long wave scheme has less spectral resolution than the 
pre-existing scheme of Fels and Schwarzkopf (1975), but it is not thought that this 
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simplification gives rise to serious errors. The short wave and long wave radiation 
tendencies are now calculated every hour and every three hours, respectively, which 
allows effective resolution of the diurnal cycle. 

2. The description of the aerodynamic resistances (or surface drag coefficient) for 
momentum and sensible/latent heat transfer between the lowest model layer and 
the surface was based on the Monin— Obukhov similarity theory as used in the 
E2— physics package of the Geophysical Fluid Dynamics Laboratory (GFDL) model 
(Miyakoda and Sirutis, 1986). When the virtual temperature correction is 
neglected, the surface drag coefficient depends on only two nondimensional 
parameters, the Richardson number and z r /z 0 , where z r is a reference height located 
within the lowest model layer and z 0 is the surface roughness length. A simple 
analytic function was developed that gives results dose to those calculated with the 
original formulation (Appendix A). This function is used in order to avoid the 
costly iterations necessary for calculating the Monin-Obukhov length. 

3. The level 2 second-order closure model of Mellor and Yamada, (1982) was 
implemented as a parameterization of vertical diffusion in the PBL. Tins has been 
used in the GFDL spectral model (Miyakoda and Sirutis, 1977) and the fine mesh 
grid model of the Japanese Meteorological Agency (Yamagishi, 1980). In this 
version of the scheme, the eddy transfer coefficient is a diagnostic function of 
Richardson number which makes it computationally expensive in operation. 

4. The surface and PBL physics were coupled so that heat and mass exchanges 
between the surface and atmosphere satisfy energy and mass conservation laws (see 
Sec. 4). An implicit scheme with explicit coefficients was used with iteration to 
time-integrate the coupled system. 

This modified version of the NMC GCM was used as the control model (Ctl— GCM) 
in this study. 

The hydrological model in place in Ctl-GCM is the same as that used in the NMC 
and GFDL GCM’s. The scheme consists of a conceptual bucket for each land point which 
is filled by precipitation and emptied by evaporation and runoff. The rate of evaporation is 
determined by the product of a potential evaporation rate; that is, the rate that would 
result from a saturated surface with the same roughness and surface temperature; and a 
wetness value, which is taken as a simple linear function of the level of water in the bucket. 
The scheme has the attraction of being conceptually simple, but cannot be considered a 
realistic description of the energy partition process as it occurs in nature. The most 
important omission in the bucket model is that of biophysical control of transpiration: 
vegetation can exert a considerable resistance on the transfer of water from soil to 
atmosphere and this can reduce the evapotranspiration rate to a level significantly below 
that calculated with the bucket formulation. Sellers (1987) discusses this and related 
shortcomings of the bucket model and concludes that it will almost always overestimate 
evaporation rates over the land. In this paper, we investigate the effect of replacing the 
conventional bucket hydrology model used in Ctl— GCM with the Simple Biosphere model 
(SiB) of Sellers et al. (1986). 

SiB is described in Sellers et al. (1986) and is reviewed briefly in Sec. 3 of this 
paper; basically, it is a model which directly addresses the effect of vegetation on land 
surface— atmosphere interactions, by modeling those physiological and physical 
(biophysical) processes which influence radiation, momentum and heat transfer. This may 
be compared with the abiotic schemes currently used in most GCM’s, see Sellers (1987). 
Here, we investigate the effect of replacing the conventional bucket hydrology model used 
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in Ctl— GCM with the Simple Biosphere model (SiB) of Sellers et al. (1986). SiB directly 
addresses the effect of vegetation on land surface— atmosphere interactions by modeling 
those physiological and biophysical processes which influence radiation, momentum, mass 
and heat transfer 

Implementation of SiB necessitated a number of modifications to the control version 
of the GCM, as discussed in Section 3 of this paper. The complete Ctl-GCM combination 
was then used to produce 30-day summer and winter simulations which are compared to 
equivalent simulations produced by the Ctl— GCM, see Sec. 4. In both SiB— GCM and 
Ctl-GCM, the same boundary conditions for sea surface temperature (SST), sea ice 
distribution, initial snow amount and topography were used. The other land surface 
boundary conditions or parameters were treated differently in the two models. 

Surface albedo in Ctl-GCM is interpolated from the seasonal climatology of Posey 
and Clapp (1964). In SiB-GCM, the albedo is calculated with a vegetation canopy 
radiation model based on the two stream approximation (Sellers, 1985). The albedo 
depends on vegetation type and state and the solar zenith angle (see Dorman and Sellers, 
1989). Figure 1 shows the difference between the SiB-GCM and Ctl-GCM albedo fields 
for the June— July period. Over most of the continental areas except for north of 50* N the 
Posey and Clapp (1964) albedo field is unrealistically low for reasons discussed in Dorman 
and Sellers (1989). The high albedo to the north of 65* N in Ctl-GCM is appropriate for 
snow— covered conditions. 



0 60E 120E 180 120W 60W C 


Figure 1: Difference map of the global albedo fields as generated by SiB-GCM and as 

prescribed for Ctl-GCM from the data of Posey and Clapp (1964). The 
fields are averages for the June-July Period. The solid contours indicate 
regions where SiB-GCM generates higher albedo values, the broken contours 
indicate where Ctl-GCM has higher albedo values. In all cases, the contour 
interval is 0.05. 
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Surface roughness length, z 0 , is time invariant in Ctl— GCM while in SiB— GCM it 
depends on the vegetation type and month of the year (see Dorman and Sellers, 1989). The 
major differences between the two models are to be found in the deserts of China and 
Mongolia: the roughness lengths there are around 30 cm in Ctl— GCM and 1 cm in 

SiB— GCM which may be compared to values of 1 cm reported by Matthews (1984). 

Up to now, the normal practice in Ctl-GCM has been to update the global soil 
moisture field at daily intervals during a simulation, thus effectively maintaining the 
wetness at or near its ‘climatological’ value at every grid point. This procedure was 
retained for the Ctl— GCM runs discussed in this paper, but a test was executed in which 
Ctl— GCM calculated the time— evolution of the soil moisture. It was found that there were 
only small differences between the wetnesses and heat flux fields produced by those two 
versions of the Ctl-GCM (see later sections). In SiB-GCM, soil moisture was predicted 
throughout all simulations. 

The soil moisture and temperature initializations are reviewed in Appendices D and 
E respectively. As discussed in Sellers et al. (1986), the initialization of the soil moisture 
field is particularly problematical as it has a relatively long ‘memory’; the soil column in 
mid latitude zones is capable of storing several hundred mm of water which translates to 
one or two months worth of evapotranspiration and runoff. In Appendix D, we describe 
the technique for initializing the SiB-GCM soil moisture field so as to make it equivalent 
to the Ctl-GCM initial field; equivalent, that is, in that a given SiB-GCM grid point is 
assumed to have been exposed to more-or— less the same climatic history as its Ctl-GCM 
counterpart up until the start of the simulation run. This technique allows us to transform 
the Ctl-GCM initial soil moisture fields into equivalent SiB-GCM fields. (The former 
were prepared by Willmott et al. (1985) using an hydrological budget methoa forced by 
observations of precipitation and air temperature.) The accuracy of these initial fields is 
questionable, a fact which is all the more disquieting as (a) it is certain that soil moisture 
initialization has a profound effect on the ensuing time-series of surface heat flux fields and 
(bl the simulation runs discussed in this paper are comparatively short. The obvious 
solution to the problem is to run both GCM’s for several years to generate self-consistent 
soil moisture climatologies, an option that was not available to us during this study due to 
constraints on computer time. 

A final comment on the effects of the soil moisture initialization: within the upper 
range of soil wetness values; that is, greater than about 0.5, both models’ responses are 
almost completely insensitive to the soil wetness, see Sellers and Dorman (1987). It is only 
when the soil moisture content drops below these values that we expect an impact on the 
calculated evapotranspiration rate. It will be shown in later discussion that the greatest 
differences between the land surface evaporation rates calculated by the two models 
occurred in zones with relatively high soil moisture contents. 

The two GCMs, Ctl-GCM and SiB-GCM, were integrated for 50 days from 
00 UTC June 15, 1986, and 30 days from 00 UTC December 15, 1985. The simulated fields 
were compared to investigate the impact of including a biophysically— based surface 
parameterization in a GCM. Particular attention has been paid to the differences between 
the land surface hydrology and energy budgets produced by the two models. 

Two additional 30— day summer case experiments were made: A SiB-GCM run 
with a reduced soil moisture initialization in North America and a Ctl-GCM run in which 
the time-evolution of the soil moisture was calculated rather than prescribed from 
climatology. 
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Following the completion of all the experiments, four coding errors were found; one 
in Ctl— GCM, one in SiB— GCM and two common to both models. In Ctl— GCM, the surface 
sensible heat flux was estimated by integrating the time rate of change in virtual 
temperature (rather than actual temperature) over the atmospheric column. This 
diagnostic error leads to errors in the Ctl-GCM sensible heat flux fields on the order of 1%. 
In SiB— GCM, the moisture contribution to the time— tendency of the virtual temperature 
due to vertical diffusion was treated incorrectly which again is thought to lead to errors on 
the order of 1%. Third, eddy diffusivities in SiB— GCM were four times as large as 
equivalent values obtained with the original Mellor-Yamada (1982) code. Last, in 
implementing the Mellor and Yamada (1982) scheme, the interlayer Richardson numbers 
were multiplied by a factor of two (relative to the original scheme) which leads to more 
rapid transitions from stable to unstable transfer conditions within the atmospheric 
column. This last error is common to both GCMs, but is not thought to have a severe 
impact on the simulations as for most of the time the transfer coefficients are operating on 
the ‘plateau’ regions (i.e., stable or unstable conditions) where the eddy coefficients are not 
sensitive to the value of the Richardson number. 

3. OVERVIEW OF THE SIMPLE BIOSPHERE MODEL (SiB) 

Here we briefly describe the philosophy, structure, atmospheric boundary 
conditions and prognostic variables of SiB. This is followed by a description of the 
implementation of SiB into the GCM and a discussion of the anticipated impact of SiB on 
the calculation of the surface heat, moisture and energy budgets. 

3.1 Model Philosophy 

In most GCMs, including the original NMC GCM and Ctl— GCM, the fluxes of 
radiation, heat (sensible and latent) and momentum across the land surface-atmosphere 
boundary are conceptualized as largely separate processes, so that a given grid point has an 
independently specified albedo, an independently specified roughness length and an 
independently formulated dependence of evapotranspiration on soil moisture. This last 
generally takes the form of the bucket model described in the preceding section. 

The interactions between the land surface and the atmosphere which are of direct 
climatological significance on time scales of seconds to years are discussed in Sellers 
et al. (1986). They may be briefly summarized as follows: 

(i) Radiation Absorption: The spectral properties of leaves and multiple 

reflections between them make vegetative canopies highly absorbent in the 
visible (photosynthetically active radiation — PAR) wavelength interval 
(0.4-0.72 fj,m) and moderately reflective in the near-infrared region 
(0.72 — 4.0 /un). In contrast, bare ground generally exhibits a relatively flat 
variation in reflectivity over the interval of 0.4 — 4.0 fa n. 

(ii) Biophysical Control of Evapotranspiration: Higher plants regulate the 

amount of gas exchange (and hence water loss) between their saturated 
interior tissues and the external air by means of valve— like structures on the 
leaf surface (stomates). These appear to respond to the environmental 
conditions of PAR flux density, temperature, humidity and leaf water 
potential so as to maximize the ratio of CO 2 influx to water vapor flux (see 
Farquhar and Sharkey, 1982) or to conserve moisture in times of water 
stress. Besides the active physiological control of evapotranspiration, plants 
also influence the rate at which free water on the land surface is returned to 
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the atmosphere: vegetative canopies intercept precipitation which can be 

directly re-evaporated to produce relatively high but short-lived releases of 
latent heat into the atmosphere. 

(iii) Momentum Transfer: Vegetative canopies usually present a rough, porous 
surface to the near— surface airflow. The resultant turbulence enhances the 
transport of sensible and latent heat away from the surface while exerting a 
drag force which may be significantly larger than that produced by bare 
ground. 

(iv) Soil Moisture Availability: The depth and density of the vegetation root 
systems and the physical properties of the soil determine the amount of soil 
moisture available for evapotranspiration. 

(v) ‘Insulation’: The soil surface under a dense vegetation canopy intercepts less 
radiation and may also be aerodynamically sheltered. As a result, the energy 
available to the covered soil is small and the component terms of the soil 
energy budget (evaporation, sensible heat flux and ground heat flux) are 
correspondingly reduced. 

In designing the SiB model, the philosophy was to model the vegetation itself and 
thereby calculate the radiation, momentum, heat and mass transfer properties of the 
surface in a consistent way. The morphological and physiological characteristics of the 
vegetation community at a grid point are used to derive the coefficients which govern the 
fluxes between the surface and the atmosphere. All of these fluxes depend upon the state 
of the vegetated surface and the atmospheric boundary conditions. 

In SiB, the world’s vegetation communities are classified into twelve types or 
biomes, see Fig. 2. Detailed descriptions of the methods used to transform the 
morphological and physiological characteristics of a given vegetation community into the 
model parameters are given in Sellers et al. (1986) and the data sets and methodologies 
used to define the global distribution of vegetation types and attributes in SiB— GCM are 
reviewed in Dorman and Sellers (1989). 

3.2 Model Structure. 

In SiB, vegetation is divided into two morphological groups: Trees or shrubs which 
constitute the upper story or canopy vegetation, and the ground cover which consists of 
grasses and other herbaceous plants. Either, both or neither of these vegetation covers may 
be present in a given grid area. The upper story vegetation consists of perennial plants 
with persistent roots assigned to a fixed depth taken to be the bottom of the second soil 
layer. The ground cover is made up of annual plants and may have a time— varying root 
depth. There is a upper, thin soil layer (soil layer 1), from which there can be a significant 
rate of withdrawal of water by direct evaporation into the air when the pores of the soil are 
at or near saturation. Beneath the root zone (soil layer 2), there is an underlying recharge 
layer (soil layer 3) where the transfer of water is governed only by gravitational drainage 
and hydraulic diffusion. The parameters required for each vegetation type in SiB are listed 
in Table 1 of Sellers et al. (1986). Values for many of the parameters are given in Dorman 
and Sellers (1989). 
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Distribution of SiB vegetation types on the 128 x 102 Gaussian grid used in the 
NMC GCM. 


1— broadleaf evergreen 

2— broadleaf deciduous 

3— mixed broadleaf 
deciduous and 
needleleaf evergreen 

4 — needleleaf evergreen 

5— needleleaf deciduous 

6— grass and broadleaf 
deciduous shrubs 

7— grass 


8 — broadleaf deciduous, shrub, 
growth singly, in patches or 
groups 

9— broadleaf deciduous, shrub, 
growth singly, or in patches 
(semi-desert) 

T— tundra 
D-desert 

C— cultivated land represented by 
wheat 
*— glacier 


The sites used for the diurnal variation studies in Sec. 4.1 are marked with 
solid dots. Energy and moisture budget analyses were made for the four 
areas encircled by solid lines, see Sec. 4.2. 





3.3 Atmospheric Boundary Conditions for SiB 

The upper boundary conditions by SiB are as follows: 

Air temperature, vapor pressure and wind speed of the lowest model layer; 

Tr> c r , Up- 

Solar zenith angle, (i. 

Five components of the incident radiation; F(o) 

F s,b(o)’ 

F s,d(o)’ 

F n,b(o)’ 

F n,d(o)’ 

F t,d(o)’ • „ „ 

Large scale and convective precipitation rates; Pl, Pc. 

3.4 Prognostic physical state variables of SiB and their governing equations 

SiB has eight prognostic physical— state variables: three temperatures (one for the 
canopy vegetation, T c , one for both the ground cover and the soil surface, T gs , and a deep 
soil temperature, Td); two interception water stores (one for the canopy, M c , and one for 
the ground cover, M g ); and three soil moisture stores (Wi, W2, and W3). The subscripts 
c and g refer to canopy and ground cover vegetation, respectively, and s refers to the 
soil surface. The combined subscript gs refers to the ground (ground cover vegetation and 
soil) and d refers to the deep soil. The subscripts 1, 2 , and 3 refer to soil layers. 

The governing equations for the three temperatures, T c , Tgs, and Td are 


Visible or PAR (< 0.72 /an) direct beam radiation, 
Visible or PAR (< 0.72 /on) diffuse radiation, 

Near Infrared (0.72 — 4.0 /cm) direct beam radiation, 
Near Infrared (0.72 — 4.0 /an) diffuse radiation, 
Thermal Infrared (8.0 — 12.0 /xm) diffuse radiation, 


Cc Tc 

— Rn c He — AEc, 

(1) 

A T gs 

= lings - H gs - AE gs - M (T gs - Tj), 

(2) 

A T d 

__ 2(Riigs — Hgs AEgs) 

— > 

(3) 





where 


Tc, T gs , Td 
Rn c , Rn g s 
H c , H gs 

He, Hgs 

Cc 

A 

r 

A 


temperature, K 
absorbed net radiation, W/m 2 
sensible heat flux, W/m 2 
evapotranspiration rate, kg/m 2 /s 
canopy heat capacity, J/m 2 /K 
latent heat of vaporization, J/kg 
day length, s 

effective heat capacity of soil, J/m 2 /K. 


The dot superscript ( • ) represents partial differentiation with respect to time. 
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In Eqs. (2) and (3), the original ‘slab’ model used to describe soil heat flux in SiB 
has been replaced by the force— restore formulation of Deardorff (1977). E c , the evapo— 
transpiration from the canopy, has two components: i) E WC) evaporation of water from the 
wetted fraction of the canopy (canopy interception loss), and ii) Ed c , transpiration of soil 
water extracted by the root system and lost from the dry fraction of the canopy. 

E gs , the evapotranspiration from the ground cover and the surface soil layer has 
three components: i) E wg and ii) Edg which correspond to E wc and Ed C for the ground 
cover, and iii) E s , direct evaporation of water from the surface soil layer. 

The governing equations for the two interception water stores are 


M c 

— Pc Dc E w c/Pwj 

(4) 

and 

M g 

= Pg — Dg — E wg / p w 

(5) 

where 

M c , M g 

= water stored on the leaves, m 


Pc Pg 

= rate of precipitation interception, m/s 


D C) Dg 
Ewc, E wg 

= water drainage rate, m/s 

= rate of evaporation of water from the wet 

portions of the 

leaves, kg/m 2 /s 

Pw 

= density of water, kg/m 3 



The governing equations for the three soil moisture stores are 


W, 

and 

W 2 

and 

W 3 


where 

W h W 2 , W 3 

6s 

Di 
Q i»i +1 

q 3 

Edc>i) Edg,i 

Pi 


K s 

Po 


p 


= ^ [ p i - Qi» 2 - ^ (E s + Edci + E dg ,i)] (6) 

= ^U^[Q>»a“Q2,3-^(Edc,2 + E dg , 2 )] (7) 

= [$2,3 - Q 3 ] (8) 

= soil moisture wetness in the three soil layers 
= 6i/Ss 

= volumetric soil moisture in layer i, m 3 /m 3 
= value of 6\ at saturation, m 3 /m 3 
= thickness of ith soil layer, m 
= flow between i and i + i soil layer, m/s 

= gravitational drainage from recharge soil moisture store, m/s 
= canopy, ground cover abstraction of soil moisture by transpiration 
from the ith soil layer, kg/m 2 /s 

= infiltration of precipitation into the upper soil moisture store, m/s 
= minimum (P 0 , K s ) 

= 0, when Wi = 1 

= hydraulic conductivity of saturated soil, m/s 
= effective precipitation rate on soil surface, m/s 
= P - (Pc + P g ) + (D c + D.) 

= Precipitation rate above the canopy, m/s 
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The precipitation excess, (P 0 — Pi), joins the gravitational outflow from the lowest 
soil moisture store, Q3, to produce runoff, Q r > 

Qr = P0-P1 + Q3 ( 9 ) 

Equations (1) through (8) represent the governing equations for SiB and may be 
used in conjunction with time series of meteorological data to simulate the energy balance 
of a given site. 

3.5. Implementation of SiB into the NMC Spectral Model 

Our experience to date has shown that it is necessary to adapt the implementation 
of SiB into the individual GCM in question. As first implemented in the Goddard Labora- 
tory for Atmospheres (GLA) Fourth-order GCM, all of the prognostic variables except for 
T c and Tg S were time-integrated with a forward differencing scheme. The two ‘fast’ 
variables T c and Tgs (which have relatively rapid time rates of change due to the large 
ratio of the flux terms fRn, H, AE) to the heat conduction/storage terms (scaled by C c 
and A) were integrated Dy a backward implicit scheme which assumes that T r , e r and u r 
are constant through the time step. This scheme produced no numerical ill-effects as in 
the GLA GCM, the lowest model layer is very thick (-100 mb) which makes T r and e r 
relatively insensitive to the addition of sensible and latent heat fluxes over short (< 1800 
secs) time steps. The backward differencing technique used to solve (1) and (2) therefore 
involved taking partial derivatives of the right hand side of (1) and (2) with respect to T c 
and T gs and solving for AT C and AT gs , the net changes in T c and T gs over a time 
step. 

Implementation of SiB into the University of California at Los Angeles/ Goddard 
Laboratory for Atmospheres (UCLA/GLA) model of Randall et al. (1985) could not be 
achieved so simply. The UCLA/GLA model has a planetary boundary layer (PBL) 
formulation in which the PBL height, 7 tb, may vary from a few millibars (~ 10 m) to over 
180 mb (- 1.5 km). The time rates of change of the PBL potential temperature, 0 b , and 
the PBL mixing ratio, are highly variable, being large when the PBL is shallow, and small 
when the PBL is deep. The formulations governing 7tb, 0 b , and qB are complex in 
themselves as they involve surface heating terms, internal PBL dynamics and entrainment 
terms. 

A scheme was developed which ensures the solution of (1) and (2) with equations 
for the time rate of change of the PBL height, potential temperature and mixing ratio 
without numerical problems. Basically the variables 0 b and qB are added to the ‘fast’ 
variable set and appropriate time derivatives of Tgs, T c , 0 b and qB are calculated at 
the beginning of each time step, (xb is advanced using a forward differencing scheme). 
These equations are then solved using a simple backward differencing scheme. This 
implementation (in the UCLA/GLA model) was used as the basis for the SiB implemen- 
tation in the modified version of the NMC model. 

In the NMC model, the Mellor and Yamada (1982) vertical diffusion model (level 2) 
was incorporated. In this parameterization, there are no explicit prognostic variables that 
describe the state of the PBL: the whole atmosphere, including the PBL, is discretized 
into layers and the state variables such as temperature, specific humidity and wind velocity 
are defined and predicted on each level; specifically, temperature and wind are defined at 
all of the model’s 18 levels while specific humidity is defined at the 12 lowest levels. These 
atmospheric variables were coupled with the SiB prognostic variables, T c and T gs , so 
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that energy and mass are conserved and an implicit scheme was used to integrate the 
system stably in time. The other (slowly— varying) SiB variables were integrated with a 
forward scheme. 

The tendencies for the time rate of change of the temperature, 6 1 , and specific 
humidity, qi, of the lowest atmospheric layer resulting from the transfer of heat fluxes 
from the surface are 


and 

= 

F ' + C^Ipl 

(10) 

M% 

= 

F ' + ^- 

(11) 

Ml 

~ 

pT 



= 

f1 


c P 

g 

H 

= 

specific heat of air, J/kg/K 
acceleration due to gravity, m/s 2 
total surface sensible heat flux, W/m 2 


E 

Api 


total surface evaporation rate, kg/m 2 /s 
lowest atmospheric layer thickness, mb 



In order to integrate the coupled system stably in time, the equations were made implicit 
with explicit coefficients. 



Rn c — H c — AE C (14) 


f A l I 5 AE es 

l A + ?lt + OT^-J 


T gs + 


f _j_ dAE gs 



w. 


To + ^«. + ff-q,= 


Hag, - H gs - E gs -il^Tg^Td) 


( 15 ) 


12 


In Eqs. (14) and (15) no time derivatives of the surface net radiation terms, Rn c 
and Rn gS) were calculated as these are small compared to the other terms. Flux 
& d H 

derivatives like were evaluated at time level t-At and once the system was 

integrated, T c , T gs , 0 U qi at time level t+At were used to compute the new values of 
the non-neutral aerodynamic transfer coefficients. The four equations were then used to 
integrate Tc> Tgs, 0 U qi again. The time tendencies for the first layer were obtained by 

dividing — 0^^ and q t+i ^ t — q t— by 2A t . Time tendencies of potential 

temperature and specific humidity at other levels were computed by solving the matrices of 
the diffusion equations successively from the bottom to the top of the atmosphere 
(Appendix B). These tendencies were added to the tendencies due to advection and other 
physical processes to give the total tendency for each prognostic variable. 

In treating the vertical diffusion of momentum, an implicit noniterative time 
integration scheme with explicit coefficients was used. The surface drag coefficient was 
evaluated using T c , T gs , 0\, qi at time level t+At after the first integration of Eqs. (12) 
to (15). 

Ctl-GCM uses the same basic time integration scheme. In Ctl — GCM, the 
atmospheric variables (temperature at 18 levels and specific humidity at 12 levels) were 
coupled with the ground surface temperature and soil temperatures at 10 cm and 50 cm 
depth. The system was integrated using the same implicit time integration scheme with 
explicit transfer coefficients. Iterations were made using a modified Newton— Raphson 
method to insure against numerical oscillations. 

As currently implemented, SiB has a significant impact on the computer time 
required for a simulation; SiB— GCM is typically 20 percent more expensive 
(computationally) than Ctl— GCM. However, it should be noted that as yet no great effort 
has been dedicated to making the code more efficient. 

3.6 Anticipated impact of implementing SiB 

The impact of implementing SiB into the GCM can be anticipated by comparing 
the surface parameterizations of SiB— GCM and Ctl— GCM. 

First, there are some differences between the albedo field calculated by SiB-GCM 
and that prescribed in Ctl— GCM, see Fig. 1. Generally, the SiB-GCM albedos are 
significantly higher over vegetated regions, leading to reductions in the surface net 
radiation of around 5% in SiB-GCM relative to Ctl— GCM. 

Second, there are some differences between the two roughness length fields, but 
these are not thought to contribute greatly to the differences between the performances of 
the two models except for the extreme cases referred to in a previous section (China, 
Mongolia). This is because the aerodynamic resistance is inversely related to the logarithm 
of the roughness length and so the turbulent heat flux calculations are relatively insensitive 
to fairly large changes in the absolute value of this parameter, see Sellers and Dorman 
(1987). 

Third and most importantly, Ctl-GCM and SiB-GCM calculate the partition of 
energy at the surface in very different ways. The effect of this can be simply illustrated by 
comparing the energy budgets calculated by the two surface parameterizations for a 
densely vegetated area. If time— invariant forcing (T r , e r , u r , F(o)) is assumed and the 
radiation budgets for the two models are taken as approximately equal, heat storage terms 
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may be neglected and equilibrium surface fluxes and surface temperatures may be 
calculated with the two parameterizations. Under these conditions, the total energy 
budgets may be written as 


Rn 

R n 


h b + ae b 
H s + * E S 


(16a) 

(16b) 


Rn 

H 

E 

A 


= net radiation, W/m 2 
= sensible heat flux, W/m 2 

3 , kg/lK 


evaporation rate, kg/m 2 /s 
latent heat of vaporization, J/kg 

The subscripts b and s refer to the ‘bucket’ and SiB parameterizations, respectively. 

The sensible heat fluxes calculated by the two models are different due to the 
different equilibrium surface temperatures. 


H = 


H„ = 


S,B 


(T B • Tr) BO 
r a pC? 

(17a) 

(T. - T r ) 
r a /?Cp 

(17b) 


P,Cp 

r a 


C„ = 


surface temperature for SiB, bucket parameterizations, respectively 
K. 

density, specific heat of air, respectively; kg/m 3 , J/kg/K. 

aerodynamic resistance between surface and lowest model layer, 
s/m 

!/(€„ Ur) 

heat transfer coefficient. 


For Ctl— GCM, the evaporation rate is given by the so-called ‘Beta— function’. 


AE„ = P 


(18a) 


e*(T) = saturated vapor pressure at temperature T, kPa 

7 = psychrometric constant, kPa/K. 

= 0.066 kPa/K at 15° C. 

P = o.rs^Wax' ’ °-^- 1 

W = soil moisture content, mm 
W ra ax = Maximum value of W, 

= 150 mm 
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In SiB, the calculation of evapotranspiration is fairly complex, involving treatments 
for two layers of vegetation and the evaporative contribution from the bare soil. However, 
for a dense vegetation canopy, like a tropical forest, the evapotranspiration rate can be 
approximated by 



’e*(T g ) - e r 
. r a + r c . 


7 


(18b) 


r c = bulk stomata! (canopy) resistance, s/m 


The essential term in Eq. (18b) is r c , the bulk stomatal or canopy resistance of the 
dry vegetation. 

Equations (16), (17), and (18) can be combined and the surface temperatures, T 

u 

and T B , eliminated using the linear approximation for saturated vapor pressure as a 
function of temperature as shown in Monteith (1973) to yield 


AE b 

_ A Rn + 

pCu <5e/r a 

(19a) 

A + t/p 


AE S 

_ AR n + pC D £e/ra 

(19b) 

A + 7 

r a + r c 



r a J 



& 

= vapor pressure deficit, kP a 
= e*(T r ) - e r 



A = slope of saturated vapor pressure versus temperature curve at 
temperature T r , kP a /K 
= 0.11 kPa/K at 15“ C 

Equation (19b) is simply the Penman— Monteith evapotranspiration equation, see 
Monteith (1973). For the same meteorological (T r , e r , u r ) and radiative (R n ) forcings, the 
ratio of the calculated evapotranspiration rates is then given by 


E 


B 



A + 7 


r a + r c 
r a 


A + 'i/P 


( 20 ) 


For a tropical forest site, r a is typically 10 s/m and, if soil moisture is not 
limiting, (3 = 1 and r c ~ 100 s/m, see Shuttleworth et al (1984, a,b), which yields a 
value of E B /E S = 4.75! Furthermore, to obtain E B /E g = 1, @ must be reduced to 0.091 

which implies a soil wetness of W = 0.07, clearly below the range of reasonable ‘tuning’ for 
a humid region. 

The above analysis is extreme as it assumes that these is no feedback from the 
surface evaporation rate on the atmospheric forcing. In Ctl — GCM, the huge latent heat 
fluxes calculated in the first few timesteps should bring about an unrealistically cold, stable 
and humid PBL which then reduces subsequent evaporation rates to less extreme values. 
Nonetheless, we should expect to see marked differences between the surface heat fluxes 
calculated by the two models with associated differences in their near surface climates. 
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4. COMPARISON OF SIMULATIONS AS PRODUCED BY CONTROL AND SiB 

VERSIONS OF THE NMC GCM. 

In this section, the simulations generated by the SiB— GCM and Ctl — GCM runs are 
compared with each other and with observations over a range of spatial scales. 

First, analyses of surface and near— surface processes for representative grid points 
within each model (see Fig. 2) were compared with equivalent field measurements. In 
particular, the 30— aay means of the simulated hourly net radiation, sensible heat and 
latent heat fluxes for a tropical forest, a coniferous forest and an agricultural site are 
compared with observations (see Sec. 4.1). 

Second, the simulated time-series of the energy and moisture budgets of four large 
regions: the Amazon Basin, central and eastern United States, southern Asia and the 

Sahara (see Fig. 2) — are discussed with reference to observations where available 
(Sec. 4.2). 

Third, the global fields of latent heat flux (evapotranspiration), sensible heat flux 
and precipitation generated by SiB— GCM and Ctl— GCM are compared (Sec. 5.3). 

Fourth and last, the impacts of the two land surface parameterization schemes on 
the large scale atmospheric circulation are discussed (Sec. 4.4). 

4.1 Diurnal variation of surface and near— surface processes. 

Time series of fluxes and other variables for several grid points over the globe were 
examined to see how the implementation of SiB affects the diurnal variation of surface and 
PBL processes. The first 30 days of the SiB-GCM and Ctl-GCM runs were analyzed for 
this study; 30 rather than 50 days were used as the surface energy budget changed rapidly 
and in a spatially heterogeneous way in the last 20 days of the Ctl— GCM run, see later 
sections. 

4.1.1 Tropical rain forest, Amazon Basin 

Figures 3 through 9 show 30-day means of the simulated diurnal variation of 
various physical parameters at a tropical forest grid point near Manaus in the Amazon 
basin (59° W, 3°S) for the period June 15— July 14. Some equivalent observations of 
Shuttleworth et al. (1984a, b, and private comm.) are also presented for comparison. 

The difference between the albedos used by the two models may be seen in Fig. 3. 
In SiB-GCM the albedo is computed using a two stream approximation model (Sellers, 
1985) which calculates the dependency of the albedo on the solar zenith angle. The 
minimum albedo at noon, (11.5%), and the daily mean albedo (12.6%) as calculated by 
SiB-GCM agree well with the observations made over the Amazon forest (Shuttleworth, 
1984b); 11.0% and 12.5% respectively. The albedo used in Ctl— GCM (Posey and Clapp, 
1964) is a climatological one and has a seasonal variation. It is unrealistically small for 
reasons discussed in Dorman and Sellers (1989). 

Because the albedos are different, there is a slight difference between the net 
radiation values calculated by SiB-GCM and Ctl-GCM (Figs. 4a, b); but the peak 
observed net radiation is 50-100 W/m 2 less than both of the simulated values. However, 
the observed net radiation is always significantly smaller than the sum of the observed 
latent and sensible heat fluxes which suggests that there may be some errors in the 
observed heat fluxes and/or net radiation values. Additionally, it is likely that the fixed 
zonal mean climatological cloudiness used in both GCM’s leads to errors in the simulations. 
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Energy Flux (W/m 2 ) 



Figure 3: 30— day mean (June 15 — July 14) 

surface albedo at a tropical forest 
grid area near Manaus (59’ W, 
3° S) as simulated by SiB— GCM 
(solid line) and as given by Posey 
and Clapp (1964) for Ctl— GCM 
(dashed line). The period June 15 — 
July 14 is used for all the 30— day 
summer means referred to below. 


SiB Rn 

• Observed Rn 

Control Rn 

• 

Observed Rn 

SiB H 

+ Observed H 

Control H 

+ 

Observed H 

SiB \E 

0 Observed E 

Control \E 

o 

Observed \E 




Figure 4: 30— day surface energy balance for the tropical forest grid point near Manaus as 
simulated by SiB and as observed by Shuttleworth et al. (1984a). The lines on 
the figure represent simulation, the points (circles and crosses) represent the 
means of observations taken over the period June 15, 1984 — July 14, 1984 at a 
site 30 km from Manaus, Brazil. a) SiB— GCM b) Ctl— GCM 
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The heat flux into the soil is not shown in Figs. 4a and 4b. Typically, this term is 
calculated to be around 20 W/m 2 during the day and —20 W/m 2 at night. The daily mean 
value is very small and is an insignificant part of the surface energy budget in both 
models. 

Figures 4a and 4b dramatically illustrate the improvement of the simulated surface 
energy budget due to the implementation of SiB— GCM. The simulations may be directly 
compared with the observations: 30 days (June 15 — July 14, 1984) of the Shuttleworth 
et al. (1984a) data set were averaged to produce the observed hourly fluxes shown in 
Figs. 4a and 4b. In SiB— GCM the maximum values of latent heat flux, 380 W/m 2 , and 
sensible heat flux, 170 W/m 2 , agree fairly well with these observations. In Ctl — GCM 
(Fig. 4b), however, almost all of the net radiation is used for evaporation while the sensible 
heat flux is almost zero. 

This unrealistic partitioning of net radiation in Ctl— GCM is a direct consequence of 
using the bucket model parameterization of surface evaporation. In the bucket model, the 
atmosphere is directly connected to the surface by the aerodynamic resistance which is 
relatively small when the surface layer is neutral or unstable, see Eq. (18a). In SiB— GCM, 
the transpired water also has to traverse the stomatal resistance wnich is typically one 
order of magnitude larger than aerodynamic resistance over the Amazon forest, see 
Eq. (18b). (In both SiB— GCM and Ctl— GCM, sensible heat fluxes are only regulated by 
the aerodynamic resistances). The morning and midday values of the aerodynamic 
resistance, r a , are around 10 s/m in both SiB— GCM and Ctl— GCM, but the stomatal 
resistance, r^, plus leaf boundary layer resistance, rb, in SiB— GCM never drops below 100 
s/m. According to Shuttleworth et al. (1984a), daytime values of the aerodynamic 



Figure 5: 30— day mean (June 15 — July 14) 

surface (stomatal) and aerodynamic 
resistances at Manaus simulated by 
SiB — GCM and Ctl— GCM. r^ is 

aerodynamic resistance; r c is 
surface (stomatal) resistance for the 
forest canopy in SiB. There is no 
equivalent to r c in Ctl— GCM. 
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• Observed T m 


SiB T c 
SiBT m 



Figure 6a: 30— day (June 15 — July 14) mean 

surface and soil temperatures at 
Manaus simulated by SiB— G CM. 
The solid line is for the tempera- 
ture at the lowest model atmos- 
pheric level, T m (<7 = 0.995) the 
coarse dashed line for the canopy 
temperature, T c , and the dotted 
line is for the ground surface/cover 
temperature, Tg$. The series of 
solid dots represent the 30— day 
mean hourly air temperatures 45 m 
above ground (10 m above canopy 
top) as observed by Shuttleworth 
et al. (1984a) at Manaus for the 
period June 15, 1984 — July 14, 
1984. 


• Observed T m 

Control T s 

Control T m 



Figure 6b: 30— day (June 15 — July 14) mean 

surface and soil temperature at 
Manaus simulated by Ctl— GCM. 
The solid line is for the tempera- 
ture at the lowest model level 
(a ~ 0.995), T,n, and the dashed 
line line is for the ground surface 
temperature, T s . The series of 
solid dots represent 30— day mean 
hourly air temperatures at 45 m 
above ground (10 m above canopy 
top) as observed by Shuttleworth 
et al. (1984a) at Manaus for the 
period June 15, 1984 — July 14, 
1984. 


resistance, r a , estimated from eddy correlation data, frequently dropped to around 10 s/m 
and the mean minimum value of the canopy resistance for the tropical forest was estimated 
to be around 130 s/m from the flux observations. The SiB-GCM simulations of r a and rb 
+ r c concur with these observations. Additionally, a ‘bump’ in the diurnal course of the 
stomatal resistance, centered on 1100—1200 local time, may also be seen in Fig. 5. This 
feature corresponds to the midday depression in stomatal conductance and photosynthesis 
as observed in Oak shrub by Tenhunen et al. (1984) and is caused by the vegetations’ 
response to the elevated temperatures and raised vapor pressure deficit in the canopy air 
space at this time of day. The reproduction of this feature indicates that the SiB-GCM is 
realistically simulating the coupled system of biological and near-surface meteorological 
processes. 
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Because of stomatal resistance, SiB— GCM generates more sensible heat than 
Ctl-GCM in most vegetated areas. This is associated with higher surface temperatures 
and a larger diurnal range in SiB— GCM. Figure 6a shows the SiB canopy, ground surface 
and the air temperature at the lowest GCM level. Figure 6b shows corresponding results 
for Ctl— GCM: ground surface temperature and the air temperature of the lowest 

atmospheric level are shown. The monthly mean air temperature 44.6 m above ground 
(average canopy top height is 35 m) as observed by Shuttleworth et al. ( 1984a, b) is shown 
in both figures. The observed temperature reaches 302 K with a diurnal range of about 
6 K: the equivalent simulated values from SiB— GCM, 302 K and 6.5 K respectively, are 
comparable to the observed values and are more realistic than the values calculated by 
Ctl-GCM (298 K and 3 K). 

The large sensible heat flux generated by SiB-GCM also helps the PBL to grow 
deeper and make the transport of heat, moisture and momentum flux in the PBL more 
efficient. Figures 7a and 7b show the calculated diurnal variation of potential temperature 
for the five lowest GCM levels (<r = 0.9950, 0.9815, 0.9605, 0.9205, 0.8565) in SiB-GCM 
and Ctl-GCM, respectively. With SiB-GCM the simulated PBL grows up to level 4 (i.e., 
the PBL depth is about 80 mb) and sometimes up to level 5 (150 mb). On the other hand 
in the Ctl-GCM simulations, PBL growth is capped at the third level which means that 
the PBL depth does not exceed 20 mb in this area. Martin et al. (1988) observed the PBL 




Local Time (Hours) Local Time (Hours) 


Figure 7: Simulated 30— day (June 15 — July 14) mean hourly potential temperatures for 

the lowest Five model levels ( O — 0.9950, 0.9815, 0.9G05, 0.9205, 0.8505) at 
Manaus. a) SiB-GCM b) Ctl-GCM 
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structure over the Amazonian rain forest near Manaus for the three week period, July 15 — 
August 7, 1985, during which the PBL depth grew to about 170 mb (1600 m) on fine days 
and 80 to 90 mb (800 m) on cloudy days. On the average, the maximum PBL depth was 
about 130 mb (1200 m) which concurs with the PBL growth simulation by SiB— GCM. 

The more unstable daytime conditions associated with the SiB— GCM simulation 
give rise to shear stresses which are roughly double those generated by Ctl-GCM, see 
Fig. 8; the 30-day mean maximum daily shear stresses are 0.28 and 0.17 Pa, respectively. 
These figures may be compared with the friction velocity measurements reported in 
Shuttleworth et al. (1984a). which are equivalent to shear stress values ranging from 0.01 
Pa (nighttime, dawn, dusk) to over 0.5 Pa (midday, unstable conditions). These higher 
shear stresses are probably associated with moderately large sensible heat fluxes, consistent 
with the SiB-GCM calculation. 


Control Stress 
SiB Stress 



Figure 8: 30— day (June 15 — July 14) mean hourly surface 

wind stress at Manaus simulated by SiB-GCM 
(solid line) and Ctl— GCM (dashed line) for the 
June— July period. 
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Kousky (1980) reported a strong diurnal pattern of precipitation in north-east 
Brazil. The time of maximum precipitation differed greatly depending on the region. It is 
interesting to see if the GCM’s generate any kind of diurnal variation in precipitation. 
Figure 9 shows the simulated diurnal variation of hourly precipitation rate and the 
observations of Shuttleworth et al. (1984a). In both simulations, precipitation is mostly 
convective with a strong diurnal variation. The simulated peak is in the afternoon in both 
cases which is not inconsistent with the observations. 


• Observed Total Precipitation 

Control Largescale Precipitation 

Control Convective Precipitation 

SiB Largescale Precipitation 



Local Time (Hours) 


Figure 9: Diurnal variation of the 30— day (June 15 — July 14) mean precipitation 

simulated for Manaus by the two GCMs for the June — July period. The solid 
line is for convective precipitation in SiB— G CM, the coarse dashed line is for 
large-scale precipitation in SiB— GCM, the dotted line is for convective 
precipitation in Ctl— GCM and the fine dashed line is for large-scale 
precipitation in Ctl— GCM. The solid dots are for equivalent observations by 
Shuttleworth et al. (1984a). 
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Figure 10: 30— day (June 15 — July 14) mean surface energy balance for an evergreen 

coniferous forest grid area (98* W, 55* N) in central Canada as simulated by 
SiB— GCM. The points marked on the Figure are for equivalent quantities as 
measured by McNaughton and Black (1973) over a Douglas fir forest in British 
Columbia on July 20, 1970. a) SiB— GCM b) Ctl— GCM 


4.1.2 Coniferous Forest 

Figure 10a (10b) shows the SiB-GCM (Ctl-GCM) simulations of latent and 
sensible heat flux for a grid point in Canada (98° W, 55° N); this point is represented as 
evergreen coniferous forest in SiB-GCM. As before, the means for the first 30 days (June 
15 — July 14) of the Northern Hemisphere summer simulation are shown. Estimated fluxes 
for one day at a coniferous forest site in Western Canada calculated from observations by 
McNaughton and Black (1973) are shown in both figures. On five other days in July 
(1970) the observed net radiation, latent and sensible heat fluxes were larger but the 
Bowen ratio was always around one during the day. The simulated Bowen ratio in 
SiB-GCM is also very close to unity during the day whereas in Ctl— GCM it is almost zero. 

4.1.3 Cultivated Area, Central Europe 

Figure 11a (lib) shows the latent and sensible heat fluxes for a grid point at (17° E, 
49° N), a cultivated area in Central Europe, as simulated by SiB-GCM (Ctl— GCM). 
Observed fluxes recorded over a barley field at Ruthe in West Germany on June 20 — 22, 
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Figure 11: 30— day (June 14 — July 15) mean surface energy balance for a cultivated grid 

area (17* E, 49* N) in central Europe. The points marked on the figure are for 
equivalent quantities as measured by van der Ploeg et al. (1980) over a barley 
field on 22 June 1977. a) SiB-CCM b) Ctl-€CM 


1979 are also shown (van der Ploeg et al., 1980). As in the case of the Amazon rain forest 
and the evergreen coniferous forest in Canada, Ctl— GCM generates far too much latent 
heat flux whereas the time series of latent heat flux in SiB— GCM compares well with 
observations. Figure 12 shows the surface albedo simulated by SiB-GCM and the pre- 
scribed Posey and Clapp (1964) (Ctl— GCM) albedo for the region; observations from a 
wheat field 50 km from Ruthe, West Germany are also shown. It appears that SiB-GCM 
generates a realistic albedo while the Ctl— GCM value is too low and does not vary with 
time. 

4.1.4 Sahara Desert 

The diurnal variation of ground surface temperature at a Saharan grid point (0° E, 
23° N) is shown in Fig. 13. SiB-GCM has a higher surface albedo which contributes to the 
lower maximum ground surface temperature compared to Ctl-GCM. 

Susskind (pers. comm.) retrieved daytime (14:30 LT) and nighttime (2:30 LT) 
ground surface temperatures from satellite data. The estimated daytime and nighttime 
temperatures around (0°E, 23° N) for July, 1979 were 324 K and 296 K, respectively, which 
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Figure 12: 30-<lay (June 15 - July 14) mean 

surface albedo at the West German 
grid point as simulated by 
SiB— GCM (solid line) and 
Ctl— GCM (dashed line). The 
crosses represent observed fluxes for 
a wheat field at Volkenrude, West 
Germany as reported in van der 
Ploeg et al. (1980) for 20-21 June, 
1979. 
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Figure 13: 30-day (June 15 — July 14) means 

of the diurnal variation of the 
ground surface temperature for a 
grid point in the Sahara Desert 
(0* E, 23° N) as calculated by 
SiB— GCM and Ctl-GCM. T d 
are the SiB ground 


and 


L gs 


surface and deep soil tempera 
tures, respectively. T 50 and T s 
are the deep soil (50 cm) and 
surface temperatures for the control 
model (Ctl — GCM), respectively. 
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gives a diurnal range of 28 K. These may be compared with the equivalent simulated 
daytime ground surface temperatures of 326.5 K in Ctl— GCM and 323 K in SiB— GCM. 
The nighttime temperatures, 310 K in Ctl— GCM and 305 K in SiB-GCM, seem to be too 
high compared with observations even when uncertainties in the retrieval procedure are 
taken into account. One possible explanation is that the fixed zonal mean clouds in the 
GCMs reduces the amount of nocturnal cooling over the deserts. Excessive downward 
sensible heat flux from the atmosphere to the surface or heat conducted from the deep soil 
layer may also contribute to the higher nighttime ground surface temperatures. 

4.2 Time tendency of regional mean energy and moisture budgets. 

Area-averaged budget analyses were made for 50 days (June 15 — August 4) of 
summer simulations by both GCMs to find out how the regional energy and moisture 
budgets changed during the course of the numerical simulations. Four continental regions 
were selected: the Amazon Basin, Central and Eastern U.S.A., Asia, and the Sahara desert 
(Fig. 2). 

In SiB-GCM, the Amazon basin is uniformly covered by tropical rain forest; the 
Central and Eastern U.S. consists mainly of cultivated land, grassland and some mixed 
forest; the Asian region consists of heterogeneous land cover including tropical rain forest, 
cultivated land and desert; and the Sahara desert consists of sand with a very sparse shrub 
cover. 

4.2.1 Energy Budget 

The daily mean surface energy budget of SiB-GCM (Ctl— GCM) for the Amazon 
Basin is shown in Fig. 14a (14b) and is summarized in Table 1. The assumption of a fixed 
zonal mean climatological cloudiness and ozone distribution results in very similar surface 
insolations for the two models. This mean value, 210 W/m 2 , is close to, but a little larger 
than the July climatological value, 195 W/m 2 , as given by Budyko (1974), but these com- 
parisons should not be extended too far in view of the prescribed cloudiness used in the 
GCMs. There is some difference between the net solar radiation fluxes calculated by the 
two GCMs because of the different albedos used; 12% in SiB-GCM and 7.7% (Posey and 
Clapp, 1964) in Ctl— GCM. This albedo difference results in a 10 W/m 2 or a 5% difference 
in the net radiation. 

The striking improvement of SiB-GCM over Ctl— GCM is shown in the partition of 
net radiation into latent and sensible heat fluxes. In the Ctl— GCM simulation, almost all 
of the net radiation is used for evaporation and the sensible heat flux is negligible. In 
SiB-GCM, 84% of the net radiation goes to latent heat flux which gives a Bowen ratio of 
0.21. This difference results from the different formulations used to calculate 
evapotranspiration in the two GCMs, see Sec. 3.6. 

The latent heat fluxes, 130 W/m 2 (SiB-GCM) and 162 W/m 2 (Ctl— GCM), cor- 
respond to 4.5 mm/day (SiB-GCM) and 5.5 mm/day (Ctl— GCM) of evapotranspiration. 
Salati (1987) summarized published estimates of the mean annual evapotranspiration rate 
for the entire Amazon Basin to give a figure in the range of 3.1 to 4.0 mm/day. If we take 
the uncertainty of this estimate and seasonal variation- into consideration, the evaporation 
rate given by the SiB-GCM, 4.5 mm/day, is in better agreement with the estimates 
obtained from observations than Ctl— GCM. 
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Figure 14a: The daily mean »urface energy 

budget calculated by SiB—GCM for 
the Amazon Basin area shown in 
Fig. 2 for a 50— day summer simu- 
lation run. Each energy band 
represents: 

1— solar radiation reflected by the 
surface 

2— net loss by longwave radiation 
3a — latent heat flux by transpira- 
tion and bare soil evaporation 

3b— latent heat flux by interception 
loss 

3a+ 3b— total latent heat flux 
4 — sensible heat flux 

The top solid line represents the 
total downward solar radiation 
flux, the second dotted line 
represents the net downward solar 
radiation, the third solid line 
represents the net downward 
radiation and the bottom solid line 
represents heat transferred to soil 
by conduction. If this lowest line 
is negative, the heat transfer is 
from the soil to the air. 
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Figure 14b: The daily mean surface energy 

budget calculated by Ctl — GCM for 
the Amazon Basin for a 50— day 
summer (June 15 — July 14) 
simulation run. Each energy band 
represents: 

1— solar radiation reflected by the 
surface 

2— net loss by longwave radiation 

3— latent heat flux 

4 — sensible heat flux 

the meaning of each line is the 
same as in Fig. 14a. 
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Table 1 The mean surface energy budget as calculated by (a) Ctl— GCM, and (b) 
SiB— GCM for the four regions outlined in Fig. 1: Amazon Basin, Central 
and Eastern United States, Sahara Desert and Asia. The means are for the 
50-day (June 15 — August 4) summer simulations in all cases. The symbols 
used in the columns represent the following fluxes (W/m 2 ). 
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(a) Mean (50-day) surface energy budget calculated by Ctl-GCM. 
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The time series of interception loss for the Amazon as calculated by SiB— GCM is 
shown in Fig. 14a and in Table 1. The time-averaged value, 36 W/m 2 , represents 28% of 
the total evapotranspiration. Shuttleworth et al. (1984a) reported an equivalent value of 
30% for their measurements taken at a site outside Manaus. An annual mean value of 35% 
was reported in Salati (1987) for a small basin (25 km 2 ) 60 km north of Manaus. Salati 
(1987) also mentions the report of Franken et al. (1982) which gives a value of 33% for a 
one— year period at another site 40 km from Manaus. In general then, the interception loss 
rate calculated by SiB— GCM seems to be realistic when compared to the observations. 
This must be partly due to the formulation used to account for the hydrological effects 
associated with the non-uniform spatial distribution of convective rainfall (Appendix C). 

Figures 14c and 14d show the daily mean surface energy budget for the central and 
eastern U.S. region. The calculated downward solar radiation at the surface is around 260 
W/m 2 which is a bit less than the climatological value of 275 W/m 2 given by Budyko 

^ and Sellers (1965). The small difference in net radiation, 20 W/m 2 , comes partly 
,he higher albedo in SiB— GCM versus Ctl — GCM, 16% vs 11%, and partly from the 
increased net longwave radiation flux associated with the higher surface temperature and 
drier PBL in SiB— GCM. As in the case of the Amazon Basin, the marked difference 
between the SiB— GCM and Ctl— GCM results is to be found in the partition of the net 
radiation into sensible and latent heat fluxes: The Bowen ratio is 0.3 in SiB— GCM and 0.1 
in Ctl-GCM. 

In Ctl— GCM, there is an increase in the sensible heat flux and a decrease in the 
latent heat flux toward the end of the 50— day summer simulation which is a result of a 
progressive reduction in the prescribed climatological soil moisture content. Figure 16b 
shows how the Ctl-GCM soil moisture was decreased from 0.88 to 0.45 over the 50 days. 
In the parallel SiB— GCM run, the equivalent (prognostic! soil wetness (Fig. 16a) decreased 
very slowly, changing from 0.70 to 0.65 in 50 days. Tne main reason for this relatively 
slow change in soil moisture is that SiB has a greater soil moisture capacity than the 
bucket model: the total depth of the moisture carrying zones is typically 1490 mm, 

porosity is 0.46 and the initial soil wetness is 0.7. These values combine to give 480 mm as 
an initial soil moisture content for SiB— GCM of which some 300 mm is available for 
evapotranspiration. This is twice the maximum soil water holding capacity (150 mm) of 
the UCLA/GLA GCM of Randall et al. (1985) and the GFDL GCM as described by 
Miyakoda and Sirutis (1986). 

Figures 14e and 14f show the daily mean energy budget for the Sahara Desert. The 
total insolation in both GCM’s decreases from 295 W/m 2 to 280 W/m 2 over the course of 
the 50-day run due to the decrease in the solar inclination. The net solar radiation is 217 
W/m 2 in Ctl-GCM and 202 W/m 2 in SiB-GCM. The 15 W/m 2 difference is due to the 
difference in surface albedo which is 31% in SiB-GCM and 25% in Ctl-GCM. However, 
this difference in the net solar radiation is approximately compensated by differences in the 
net terrestrial radiation loss so the resulting net radiation is almost the same; 125 W/m 2 in 
Ctl-GCM and 119 W/m 2 in SiB-GCM. It is interesting to note that the net radiation in 
this region comprises only 42% of the incoming solar radiation as opposed to 70-80% for 
the Amazon basin and U.S.; this is due to the high surface albedo and low atmospheric 
humidity in the desert regions. Both models give very low evapotranspiration rates in this 
area so almost all of the net radiation is returned to the atmosphere as sensible heat flux. 

Figures 14g and 14h show the daily mean energy budget for the Asian region 
indicated in Fig. 2. Because this region consists mainly of desert with relatively little 
vegetated land, the conspicuous change of surface energy partition encountered in the 
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Amazon Basin and in Central and Eastern U.S.A. is not seen so clearly here, but there is 
some increase in the Bowen ratio in SiB— GCM. 

4.2.2 Hydrological Cycle 

The simulated hydrological balance of the Amazon basin is shown in Figs. 15a and 
15b. The precipitation rate is 6 mm/day in SiB— GCM and 7 mm/day in Ctl— GCM. From 
the climatological precipitation chart (Fig. F2d of Jaeger, 1976) we see that the northern 
half of the Amazon basin is wet and the southern half is dry in summer. If an area average 
is taken, we get 4—5 mm/day for the climatological rainfall rate in June— July for the entire 
Amazon Basin. Compared with this observed value, both GCM’s, especially Ctl— GCM, 
give more precipitation. The 1 mm/day decrease in the SiB— GCM versus the Ctl-GCM 
precipitation rate corresponds to the reduction in the SiB-GCM evaporation rate. The 
same situation holds in North America and Asia. In these regions, the rainfall over land 
seems to be directly affected by the regional evapotranspiration rate. 
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Figure 15: The daily mean moisture budget calculated by SiB-GCM in the Amazon 

Basin for the 50— day summer (June 15 — July 14) simulation run. The solid 
line is for precipitation, the coarse dashed line is for evaporation, the fine 
dashed line is for runoff and the dotted line is for the atmospheric horizontal 
moisture convergence. All values are in mm/day. 
a) SiB-GCM b) Ctl-GCM 
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From Figs. 15a and 15b, we see that the precipitation exceeds evaporation by 1.5 
mm/day on the average, with the difference made up by large scale horizontal moisture 
convergence. This advected moisture contributes about 25% (SiB-GCM) and 20% 
(Ctl— GCM) to the total rainfall in the region. 

The runoff rate in SiB-GCM is 2.5 mm/day on the average. This amount exceeds 
the difference between the rainfall rate and evaporation rate by 1 mm/day and implies 
that the total soil moisture content in SiB is steadily decreasing (see Fig. 16a). Oltman et 
al. (1964) estimated the annual mean outflow from the Amazon basin to be about 2.5 
mm/day from observations of the Amazon river velocity profile and stage. There is a fairly 
considerable annual variation in this quantity, but the June to August period appears to 
bracket this mean value while showing a steadily decreasing outflow rate. This decrease in 
the outflow is consistent with a fall in the average soil moisture content for the basin as 
calculated by SiB. 

In Ctl— GCM, runoff is generated by the difference between rainfall and evaporation 
when the soil is saturated. When soil moisture was predicted (see below) it decreased 
significantly, but because soil moisture in this main Ctl— GCM run was interpolated from 
monthly climatology, the soil wetness in the Amazon basin was almost always at the 
saturated value, see Fig. 15b. 
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Figure 15c: The same as Fig. 15a except for central 

and eastern U.S. 
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Figure 15d: The same as Fig. 15b except for central 

and eastern U.S. 
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The calculated moistures budget for the U.S. region are presented in Figs. 15c and 
15d. In SiB-GCM (Ctl-GCM) the calculated precipitation rate is 4.7 mm/day (6.2 
mm/day) and the evaporation rate is 4.5 mm/day (6.0 mm/day). In both cases, the 
time— mean of the horizontal convergence of moisture by tne large scale atmospheric 
circulation is very small. This is supported by the analysis of Rasmusson (1968) who 
computed the horizontal moisture convergence over the same Central and Eastern U.S. 
region for each month of the year from 5 years of radiosonde data. He estimated the 
surface evaporation rates to be 2.7 mm/day and 3.1 mm/day in June and July, 
respectively, whereas the horizontal moisture convergences were very small for the same 
period. The evapotranspiration rate obtained from the simulations exceeds the estimates of 
Rasmusson (1968) by 50% in SiB-GCM and by 100% in Ctl— GCM. (The corresponding 
evapotranspiration rate of 4.5 mm/day in SiB-GCM seems close to the values observed for 
a forest without soil moisture stress: 5 mm/day is given by Swift (1975) as the June— July 
mean for a deciduous forest in the southern Appalachians and 4 mm/day is given by 
Luxmoore (1983) as the July-September mean for a deciduous forest in the Eastern U.S.). 

From the hydrological analyses of the Amazon and North American regions, we 
note that the higher precipitation rates over land in Ctl— GCM were associated with higher 
evapotranspiration rates over the same area and not from an increase in horizontal 
moisture convergence. This leads us to suggest that the precipitation over these regions is 
in large part controlled by the transpiration there. In order to check this hypothesis, 
another SiB-GCM experiment was performed with a reduced initial soil moisture content. 

A 30-day SiB-GCM run was performed using the same initial conditions as the 
50— day summer SiB-GCM run except that in North America the initial soil wetness was 
set to 0.45 wherever it exceeded this value. Figure 15g shows the moisture budget for the 
central and eastern U.S. region as calculated for this run. Once again, the time— averaged 
horizontal moisture convergence is almost zero and the reduced evapotranspiration rate, 
down from the 4.5 mm/day of the original 50— day SiB-GCM run to less than 3.5 mm/day 
in this reduced soil moisture case, is reflected in a corresponding 1 mm/day decrease in the 
time— averaged precipitation rate. This result strongly supports the hypothesis that the 
evapotranspiration and rainfall rates are tightly coupled in this region during the summer 
months. 

Returning to the original summer run, the runoff generated for the U.S. region by 
the SiB-GCM is 1 mm/day and is not in equilibrium with the precipitation rate minus the 
evaporation rate. As in the case of the Amazon Basin the runoff is caused by the drainage 
from the recharge zone of the SiB soil moisture model. 

The moisture budget for Asia for the summer run is shown in Figs. 15e and 15f. 
The precipitation rate in SiB-GCM (Ctl— GCM) is 7.2 mm/day (8.5 mm/day), evapo— 
transpiration 3.8 mm/day (4.7 mm/day), horizontal moisture convergence 3.4 mm/day (3.8 
mm/day) and runoff 3.0 mm/day (1.1 mm/day). The decrease in precipitation 
(1.3 mm/day) in SiB-GCM compared with that in Ctl— GCM comes partly from a decrease 
in the evapotranspiration rate (0.9 mm/day) and partly from a decrease in the horizontal 
moisture convergence (0.4 mm/day). The main characteristic to note for this region is that 
horizontal convergence is calculated to be a significant contributor to the precipitation in 
both GCM’s for this region. 
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Figure 15e: The same as Fig. 15a except for Asia. 


Figure 15f: The same as Fig. 15b except for Asia. 
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Figure 15g: The daily mean moisture budget as 
calculated by SiB— GCM for 30 
days (June 15 — July 14) in the 
eastern and central U.S. where the 
initial soil wetness was limited to a 
maximum of 0.45. Other condi- 
tions and symbols are the same as 
for the 50— day summer SiB— GCM 
simulation (Fig. 15a). 
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Figure 16a: The evolution of the soil wetness 

calculated by SiB— GCM during a 
50— day summer (June 15 — 

August 4) simulation run. The 
depth— weighted average of the 
wetnesses of the surface zone, root 
zone, and recharge zone was taken. 
The solid line is for the Amazon 
basin, the dashed line is for the 
central and eastern U.S. region and 
the dotted line is for the Asia 
region; see Fig. 2. 


Asia 

N. America 

Amazon 


Ctl-GCM: 6/15-8/4 



Figure 16b: The same as Fig. 16a except for 

Ctl— GCM. Here the time— evolu- 
tion of soil wetness follows the 
prescribed climatology of Willmott 
et al. (1985). The soil wetness in 
Ctl— GCM is derived by dividing 
soil moistures by 150*0.75 mm, 
three quarters of field capacity. 
Soil wetnesses exceeding 1.0 are set 
to 1.0 in the GCM. 
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Because runoff is less than precipitation minus evaporation in SiB— GCM, the soil 
moisture in the region increases slowly from 0.45 to 0.50 during the run (Fig. 16a). At the 
same time the climatological soil moisture prescribed in Ctl— GCM increases from 0.43 to 
0.54 (Fig. 16b). 

Lastly, it will be remembered that a decision was made early in the study to 
continuously update the soil moisture content in Ctl— GCM in line with the climato- 
logical estimates of Willmott et al. (1985). An additional 30-day Ctl— GCM run was 
performed in which the soil moisture was predicted and updated by the model. Figure 16c 
shows the predicted time evolution of soil moisture which may be directly compared with 
Fig. 16b. It may be seen that the predicted soil moisture decreases more rapidly in the 
Amazon region, but more slowly in the central and eastern U.S. as compared with the 
imposed climatology. 

In Figs. 16b and 16c, the differences between the simulated and climatological 
time-series of soil moisture are due to regional differences in the precipitation— evaporation 
ratios as simulated by Ctl— GCM and calculated by Willmott et al. (1986). However, the 
difference in the wetness values associated with the prognostic and climatological soil 
moisture fields was comparatively small and so the difference between the global fields of 
evapotranspiration and sensible heat flux produced by the two Ctl— GCM runs was slight; 
at day 30, the version with the predicted soil moisture gave rise to a sensible heat flux of 
only 10 W/m 2 less over the U.S. than that obtained with the climatological soil moisture 
field. There was virtually no difference in the heat flux fields calculated for the Amazon by 
the two Ctl— GCM runs at this time. 
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Figure 16c: 30— day (June 15 — July 14) 

time — aeries of soil moisture 

predicted by the Ctl — GCM for a 
summer circulation. This soil 
moisture evolution may be 
compared directly with the 
prescribed climatology shown in 
Fig. 16b. 
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4.2.3 Soil Heat Conduction and Temperature 

Figure 17a (17b) shows the area-averaged ground surface temperature and deep soil 
temperature (soil temperature at 50 cm depth) in SiB— GCM (Ctl — GCM) for the four 
regions. In Ctl— GCM, the displayed surface temperature was derived from the daily mean 
upward longwave radiation flux at the surface; it was found that this did not differ much 
from the daily mean of the surface temperature. The deep soil temperature in SiB— GCM 
shows some initial adjustment due to uncertainties in its initialization; however, analysis of 
the surface energy budget showed that the daily mean heat conduction into or out of the 
soil was usually less than 10 W/m 2 and so was not critical to the surface energy budget. 



SiB-GCM: 6/15-8/4 



Figure 17a: The deep soil temperature, Td, at 
00 GMT and daily mean surface 
temperature, T Sl for a 50— day 
(June 14 — August 4) period 
summer simulation run by 
SiB— GCM. The daily mean 
surface temperature was derived 
from the daily mean upward 
longwave flux at the surface and so 
is a weighted mean of the canopy 
and ground temperature. The top 
pair of lines are for the Sahara 
Desert, the second pair for central 
and eastern U.S., the third pair are 
for the Amazon and the bottom 
pair are for Asia. 6 K was added 
to the temperatures for central and 
eastern U.S. for easier display. 
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Ctl-GCM: 6/15-8-4 




Figure 17b: The same as Fig. 17a, but for Ctl— GCM. The deep soil temperature, Tso, is 

at 50 cm depth. 


4.3 Global Surface Field 
4.3.1 Summer Case 

The surface latent heat flux for the first 30 days in the summer simulation by SiB— 
GCM (Ctl— GCM) is shown in Fig. 18a (18b) and the difference in Fig. 18c. Except for the 
deserts, the latent heat flux in Ctl— GCM over the continents is usually more than 150 
W/m 2 and sometimes exceeds 200 W/m 2 in parts of North America and South Asia. In 
SiB— GCM, latent heat fluxes of more than 150 W/m 2 are observed only in South Asia. 
The difference map shows that SiB-GCM gives 25—50 W/m 2 (1—2 mm/day) less evapo— 
transpiration over most vegetated land. 

We compared the latent heat flux calculated with SiB-GCM with the global map 
of climatological latent heat flux for June as estimated by Budyko (1974). Budyko’s 
estimates were found to be less than the simulated value over all land surfaces; for 
example, his estimated latent heat fluxes in the Amazon basin and central/eastern U.S. are 
80 W/m 2 and 60 W/m 2 , respectively, which are much less than the observed values of 
90—120 W/m 2 and 90 W/m 2 for the respective regions cited in the previous subsection. 
The global map of surface net radiation (Budyko, 1974) shows values which were typically 
one half of the downward solar radiation over the land surface. The simulations from both 
GCM’s suggest that this only happens over deserts where the surface albedo is high and the 
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Figure 18c: The difference between the 30— day (June 15 — July 14) mean latent heat 

fluxes of the summer simulations: SiB—GCM minus Ctl— GCM The contour 
interval is 25 W/m 2 . The land areas where Ctl— GCM produces > 25 W/m 2 of 
latent heat flux compared to SiB—GCM have been shaded. 


atmospheric humidity is low. Table 1 shows that over the Sahara desert, about 58% of the 
incoming solar radiation is lost by reflection and longwave emission and the remaining 42% 
is available for sensible and latent heat fluxes. Over densely vegetated regions like the 
Amazon basin, only about 30% is lost (see Figs. 14a and 14b). Although the net loss of 
surface energy by terrestrial radiation is dependent on the amount and height of the 
prescribed cloud in both the climatological estimation (Budyko, 1974) and the GCMs, 
Budyko’s value does seem to severely underestimate the net radiation and consequently the 
latent and sensible heat fluxes. 

Figure 19a (19b) shows the mean daily precipitation rate simulated by SiB—GCM 
(Ctl— GCM) and Fig. 19c the difference between them. For comparison, a climatological 
field for the June-July period (Jaeger, 1976) is presented in Fig. 19d. Generally, the 
simulated precipitation pattern agrees well with the climatology, but the Ctl— GCM field is 
marked by excessive precipitation in North America and northern Eurasia. The difference 
map clearly illustrates the reduction of precipitation over vegetated land in SiB—GCM. 
There is a northward shift of the monsoon rainbelt in Asia in SiB-GCM relative to 
Ctl— GCM. 

The surface sensible heat flux for the first 30 days in the summer simulation by 
SiB—GCM (Ctl— GCM) is shown in Figs. 20a and 20b with the difference shown in Fig. 20c. 
Except for the arid zones, the time— averaged sensible heat flux in Ctl— GCM is very small 
over land and even negative over the tropical forests while in SiB—GCM, the time- 
averaged sensible heat flux is more than 25 W/m 2 almost everywhere over the vegetated 
land. The negative value observed in the defference field in northern China is a result of 
the higher albedo and wetter soil in SiB-GCM, the wetter soil being due to the northward 
shift of the rainbelt in China in SiB—GCM. 
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Figure 19a: Mean daily total (convective k large scale) precipitation in the SiB-GCM 

summer simulations for the first 30— day (June 15 July 14) period. The 
contour intervals are 1.0, 2.0, 4.0 and 8.0 mm/day. Areas with rainfall rates 
greater than 4.0 mm/day have been shaded. 



Figure 19b: The same as Fig. 19a except for Ctl—GCM. 
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Figure 19c: The difference between the 30— day (June 15 — July 14) mean daily 

precipitation for the summer simulations: SiB— GCM minus Ctl— GCM. The 
contour interval is 1 mm/day. The land areas where Ctl — GCM produces > 1 
mm/day of precipitation compared to SiB — GCM have been shaded. 



Figure 19d: Mean observed daily precipitation for June— July from Jaeger (1975). 
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Figure 20c: The difference of the 30— day (June 15 — July 14) mean sensible heat fluxes for 

the summer (June 15 — July 14) simulations: SiB — GCM minus Ctl — GCM. 

The contour interval is 25 W/m 2 . Land areas where SiB — GCM produces > 25 
W /m 2 of sensible heat flux compared to Ctl— GGM have been shaded. 


Figures 21a and 21b show the daytime PBL depth as simulated by SiB-GCM 
(Ctl-GCM). The PBL depth was derived by counting the number of layers where the 
30-day mean heating rate generated by vertical diffusion exceeded 1 K/day. This 
definition was found to agree well with the conventional method of using potential 
temperature gradients. 

Except for the deserts, the PBL depth in Ctl— GCM is calculated to be about 20 mb 
over the land which is unrealistically shallow; this is a direct result of the small sensible 
heat flux generated by Ctl-GCM. In SiB-GCM the PBL depth is more than 60 mb over 
most of the vegetated land and in the Amazon basin, the simulated depth of 100 mb agrees 
well with the observations of Martin et al. (1988). In both GCM’s the PBL grew to 
300 mb at some desert grid points. This very high value is consistent with the report of 
Blake et al. (1983) that the PBL depth in the Saudi Arabian desert could reach 350 mb. 

In SiB-GCM, the daytime PBL in the SiB-GCM is more unstable than that of 
Ctl-GCM, and so momentum transport is more efficient in the former model. This results 
in a larger surface wind stress in SiB-GCM over land. Figures 22a and 22b show the 
zonally— averaged zonal surface wind stress in SiB-GCM (Ctl-GCM). The value over 
ocean is almost the same for both GCMs and agrees fairly well with that estimated from 
surface data by Hellerman and Rosenstein (1983). Over the continents, the SiB-GCM 
value is almost twice as large as the Ctl-GCM value. 
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Figure 21a: Mean maximum PBL height in the SiB— GCM summer (June 15 — July 14) 

simulation The contour interval is 30 mb. Land areas where the mean 
maximum PBL height exceeds 60 mb have been shaded. 
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4.3.2 Winter Case 

In the winter simulation, the two 30-day runs were made with the surface drag 
parameterization of Deardorff (1972) which gave a slightly larger surface drag coefficient 
under stable and neutral conditions than the formulation based on the Monin— Obukhov 
model used in the summer simulations. This different formulation tended to reduce the 
Bowen ratio slightly, but the overall difference between the two GCMs was the same as in 
the summer simulation. 

Figures 23, 24, and 25 show the differences between the latent heat flux, sensible 
heat flux and precipitation 30— day (December 15 — January 13) mean fields as calculated 
by SiB— GCM and Ctl— GCM for the 30— day winter simulation. The SiB— GCM enhances 
sensible heat flux by up to 50 W/m 2 with comparable decreases in the latent heat flux. 
The reduced sensible heat flux over Australia calculated by SiB— GCM is associated with 
its higher desert albedo values. 

The simulated precipitation patterns from both GCM’s are in good agreement with 
the climatology except in the Amazon basin. The climatological precipitation field exhibits 
one local maxima in the region; one in the northeast and another in the southeast of Brazil. 
The SiB— GCM is associated with decreased precipitation over land (through the reduction 
of evapotranspiration) but SiB did not change the basic precipitation pattern; it seems that 
other physics parameterizations in the GCM may be responsible for this problem. 
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Figure 23: Difference field of latent heat flux as calculated by SiB— GCM and Ctl— GCM 

for the 30— day winter (December 15 — January 13) simulations: SiB— GCM 

minus Ctl— GCM. The contour interval is 25 W/m 2 . Land areas where 

Ctl— GCM produces > 25 W/m 2 of latent heat flux compared to SiB— GCM 
have been shaded. 
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Figure 24: Difference field of sensible heat flux as calculated by SiB—GCM and Ctl— GCM 

for the 30— day winter (December 15 — January 13) simulations: SiB— GCM 

minus Ctl— GCM. The contour interval is 25 W/m 2 . Land areas where 
SiB—GCM produces > 25 W/m 2 of sensible heat flux compared to Ctl— GCM 
have been shaded. 



Figure 25: Difference field of mean daily precipitation rate (large-scale plus convective) as 

calculated by SiB—GCM and Ctl— GCM for the 30— day winter (December 15 — 
January 13) simulations: SiB—GCM minus Ctl— GCM. The contour interval is 
1 mm/day. Land areas where Ctl— GCM produces more than 1 mm/day of 
precipitation than SiB—GCM have been shaded. 
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4.4 Impact on the Large Scale Atmosphere Circulation 

The impact of SiB on the large scale atmospheric circulation is briefly discussed 
here. Clearly, the shortness of the simulation runs performed with SiB— G CM and 
Ctl-GCM and the uncertainties in the soil moisture initialization (see earlier sections and 
Appendix D) prevent us from drawing definitive conclusions. Longer runs, where the 
effects of the soil moisture initialization scheme can be disregarded, must be carried out to 
fully investigate the interactions between the surface energy balance and the large scale 
circulation. 

4.4.1 Summer Case 

Figures 26a and 26b show the 200 mb zonal winds for the first 30 days of the 
summer simulation by SiB— GCM (Ctl-GCM). The weakening of the tropospheric jet in 
SiB-GCM compared with Ctl-GCM is noticeable in the northern hemisphere. Wind 
speeds of more than 30 m/s are observed in the NMC analysis (Fig. 26c) and in the 
Ctl— GCM simulated field, but not in the SiB-GCM simulated field. Although not clearly 
visible in the maps, the jet stream has also shifted slightly northward in the SiB— GCM 
simulation. 



Figure 26a: Mean sonal wind at 200 mb in the SiB-GCM rammer (June 15 - July 14) 

simulation for the first 30 days. The contour interval is 10 m/s. 
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Figure 27a: Mean eddy geopotential height at 500 mb in the SiB-<2CM for the summer 

(June 15 — July 14) simulation. The contour interval is 40 m. 





Figure 27c: Eddy geopotential height at 500 mb from NMC analyses averaged for the 

period June 15 — July 14, 1986. The contour interval is 40 m. 


The 500 mb eddy geopotential height for the first 30 days of the summer simulation 
as simulated by SiB-GCM (Ctl — GCM) is shown in Fig. 27a (27b). When compared with 
the NMC analysis the simulations by both GCM’s appear to be better in the northern 
rather than the southern hemisphere (Fig. 27c). The wavy pattern over the Eurasian 
continent in the SiB— GCM simulation is flatter than in the Ctl-XSCM simulation while the 
trough to the north of the Caspian Sea is completely missed in the SiB-GCM simulation. 

4.4.2 Winter Case 

The simulation of the circulation in the northern hemisphere is very good in both 
GCM’s when compared with the analyzed field. Ctl— GCM simulates the ridge around 
Lake Baikal better than SiB-GCM, but the ridge to the south of New Zealand seems to be 
better represented in SiB-GCM. 

5. SUMMARY AND DISCUSSION 

The Simple Biosphere Model (SiB) developed by Sellers et al. (1986) was 
implemented in a modified version of the NMC spectral model. The SiB-GCM was 
integrated for 50 days in a (boreal) summer simulation and for 30 days in a (boreal) winter 
simulation. The simulated surface and atmospheric fields were compared with those 
produced by the control GCM (Ctl-GCM) which has the same physics package as 
SiB-GCM except for the ‘bucket’ hydrological model and the three soil layers with 
predicted temperatures. 

Ctl— GCM was found to give too much evaporation over the vegetated land; 
sometimes the latent heat flux was almost equivalent to the net radiation. On the other 
hand, SiB-GCM generated sensible and latent fluxes over the land which compared well 
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with observations. The improvement in the calculated Bowen ratio in SiB— GCM comes 
from the explicit treatment of transpiration control by stomatal resistance. 

There is a secondary effect due to the generally higher land— surface albedos 
calculated by SiB— GCM as compared to the prescribed fields used in Ctl— GCM. This is 
thought to alter the total available energy by about 6% in vegetated areas with a similar 
impact on the evapotranspiration rate. The stomatal resistance effect, however, accounts 
for the rest of the reduction in evapotranspiration, typically 25% or greater over the moist 
vegetated regions. This result corresponds with those in the sensitivity study of Sellers and 
Dorman (1987). 

The increase in the sensible heat flux (25 — 50 W/m 2 ) over the vegetated land 
calculated by SiB-GCM was associated with higner surface temperatures and larger and 
more realistic diurnal ranges in the surface and PBL temperatures. The daytime PBL 
depth was 60—100 mb over vegetated land in SiB-GCM compared with the unrealistic 
value of 20 mb calculated by Ctl— GCM. A deeper PBL and more active momentum 
exchange over the land led to an increase in the surface stress calculated by SiB-GCM. 

The SiB-GCM June — July precipitation rate was 1 to 2 mm/day lower than in the 
Ctl— GCM run in the Amazon basin, the African tropics, North America, the Soviet Union 
and Europe. The monsoon rainbelt in China moved northward, but overall the preci- 
pitation there also decreased. This decrease in the precipitation rate corresponded roughly 
to the decrease in the time— averaged latent heat flux, i.e., about 25—50 W/m 2 . 

SiB-GCM still generated more rain than the climatological precipitation in some 
places over land. This may be due to inadequacies in the physical parameterizations in the 
GCM or it could be that SiB-GCM still calculates evapotranspiration rates that are too 
high because of problems in the soil moisture initialization. Over the central and eastern 
United States, the simulated horizontal moisture convergence in June— July was almost 
zero in both GCMs and was in good agreement with observations, but the evaporation rate 
and hence the precipitation rate in the region was far lower in SiB-GCM. This suggests 
that the excessive precipitation calculated over land by Ctl— GCM may be the result of an 
overestimation of the evaporation rate. An experiment was performed to explore this 
hypothesis by integrating SiB-GCM with decreased initial soil moisture in central and 
eastern U.S. (The soil moisture initialization was suspect as it is an uncertain part of any 
GCM model). The results showed that the simulated horizontal moisture convergence 
remained almost zero and the smaller evapotranspiration rate over land directly 
contributed to a reduction in the precipitation rate. 

In the course of the 50— day simulation, soil wetness in SiB-GCM changed very 
slowly. This was partly due to the large soil moisture capacity in the SiB formulation and 
partly due to the small differences between evapotranspiration and precipitation. 

In Ctl— GCM soil wetness was updated every day from values interpolated from 
monthly climatology. Over vegetated land in the Northern Hemisphere the soil moisture 
decreased rapidly over the June — August period. This brought about a gradual increase in 
the sensible heat flux and a decrease in the latent heat flux which led in turn to a reduced 
precipitation rate. 

The realism of SiB-GCM and Ctl— GCM on smaller spatial scales was also 
scrutinized. The diurnal variation and the mean value of the albedo as calculated with 
SiB-GCM were found to be satisfactory when compared with observations while the 
albedos used in Ctl— GCM were generally too low. The interception loss fraction is difficult 
to verify, but at least in the Amazon basin the SiB-GCM value of 28% of the total latent 
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heat flux seems to be reasonable. Overall, SiB— GCM appears to give a more consistent 
and realistic simulation of hydrological processes and near-surface atmospheric conditions 
than Ctl — GCM. 

It should also be noted that the differences between the values of sensible heat flux, 
evapotranspiration and precipitation calculated by the two models are maximized in the 
humid vegetated regions of the summer (growing season) hemisphere. Under the prevailing 
conditions of low relative humidity and high radiation loads the evaporation rate is very 
sensitive to the value of the stomatal resistance, which is around 100 s/m in SiB— GCM and 
effectively zero in Ctl— GCM. 

A preliminary assessment of the impact of SiB-GCM on the large scale atmos- 
pheric circulation was also made. The large-scale circulation changed only slightly in the 
northern hemisphere winter simulation. This is understandable because of the small land 
area in the southern hemisphere and because the continents in the northern hemisphere 
receive little solar radiation in winter. 

In the summer simulation there was some noticeable impact: the 200 mb jet stream 
and the wind speed at 850 mb were found to be unfavorably weakened in SiB-GCM which 
also failed to simulate the stationary trough in the 500 mb geopotential height. We do not 
yet understand the causality of this effect. Clearly, the above comments must be regarded 
as very preliminary in view of the shortness of the simulation runs (50 days or less) and the 
uncertainties in the soil moisture initialization. 

In summary, the implementation of SiB in the GCM has improved the simulated 
partition of absorbed energy at the land surface. Generally, the explicit recognition of 
stomatal resistance terms in the SiB latent heat flux formulation has resulted in a 
reduction in the land surface evapotranspiration rate compared to the bucket model 
calculation. In the summer simulations, this led to approximately matching reductions in 
the continental precipitation rates in SiB-GCM. The reduction in evapotranspiration and 
rainfall rates over the land due to SiB seems to be supported by the available observations. 
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APPENDIX A 

ANALYTIC FUNCTION FOR C m AND C h BASED ON 
MONIN-OBUKHOV SIMILARITY THEORY. 


C m 


1 r _ 1 


(A- 1) 


where C m , Ch are multiplied by the magnitude of the wind velocity to yield the drag 
coefficients for momentum transfer and heat/moisture transfer. 


F° = Fg = [ln(z r /z o )/0.40] * 0.87581 - 1.5630 

stable region : Ri>0 

F m /F° = Fh/Fg = exp[— 1.8900 + (1.8900)2 + 5.0519x] 
x = ln(1.0 + 10.815 Ri) 

unstable region : Ri>0 

F m /F° = exp[ 1.2270 - (1.2270)2 + 1.2642x] 
ln(1.0 - 1.2743 Ri) 


x = 


Fh/Fg = exp[1.3462 - (1.3462)2 + 1.8380x] 

x = ln(1.0 - 3.4805 Ri) 


(A.2) 

(A.3) 


(A.4) 


(A.5) 


where F° and F£ are the neutral values of F m and Fh, z r is the reference height, z 0 is 
the roughness length and Ri is the Richardson number. 

The curves derived using the original Monin— Obukhov formulation and those using 
the analytic function above are presented in Fig. A.l. In all cases the reference height z r 
was set to 20 m or to half the height of the lowest model level. 
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APPENDIX B 

DERIVATION OF EFFECTIVE PBL MASS AND FORCING 


The vertical diffusion tendencies for potential temperature, 0 , and specific 
humidity, q , are given by 


<90 _ 1 d f v 501 

ITt ~ ~pUz l pKh ~dl] 

$< = jm 


(B.l) 


(B.2) 


where Kh is the vertical diffusion coefficient for heat and moisture transfer. The two 
equations are discretized in the vertical direction using a — coordinates. 


Tpp = Ajc (0k-i-0k)-Bk(0k-©k+i) 

77^ = Ak (qk-i— Qk)— Bk(qk— Qk+i) 


; k = 1, 2, , 18 


; k = 1,2, 


(B.3) 

(B.4) 


where Ak and Bk are coefficients which arise from the vertical discretization scheme. 
If an implicit scheme with explicit coefficients is used 


-Ak©k-i + (1 + Ak + Bk)0k - Bk©k+i = 

a A t t At -p t At pvt Atl 

Ak[0 k _i - 0 k J- B k[ 0 k - 0 k+lJ 


(B.5) 


Ok = 


Ok +At — 0k At 


A* = 2AtA , B* = 2AtB, 


Then the system is solved by elimination. 


58 


APPENDIX G 

TREATMENT OF SPATIALLY NON-UNIFORM PRECIPITATION IN SiB. 


A simple scheme is used in SiB to account for the hydrological effects resulting 
from the spatial non— uniformity of convective precipitation. In nature, the net effect of 
the ‘clumped’ precipitation patterns typical of convective storms is to increase the 
area-averaged canopy throughfall and surface (infiltration excess) runoff rates relative to 
the equivalent rates produced by a spatially uniform rainfall fields. In most GCM’s, a 
single (area— averaged) figure for convective precipitation is produced for each grid— area for 
each time step. In SiB, we assume that this rainfall amount is spatially distributed as 
shown in Fig. C.l and given by Eq. C.l. 

I c (x) = ace -bcX + c c (C.l) 

I c (x) = relative amount of convective precipitation as a function of 

fractional area of grid area x, (0<x<l), falling within a time step 

a c ,bc,c c = constants 


Figure C.l: Precipitation area— amount 

relationships used in SiB. The 
variable x refers to fraction of the 
grid area, the variable I( X ) refers 
to the relative amount of precipj — 
tat ion. Note that the large scale 
precipitation, Ip(x)» almost 

invariant over the grid area while 
convective precipitation, Ic(x)» is 
non— uniformly distributed. 


0 0.5 1 

x — ► 
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The constants a, b, c were obtained by comparison with the data of Ruprecht and Gray 
(1976) and are normalized such that 

(C.2) 

The values chosen for convective precipitation were a c = 20, b c = 20, c c = 0.206E— 8. For 
precipitation generated by large scale systems, the equivalent values are a p = 0.0001, 
b p = 20 and c p = 0.9999 which means that in the interval 0 < x < 1, I p (xJ = 1. Both 
types of precipitation (most models can produce both types simultaneously) can then be 
combined to give a single area— amount function by: 

(C.3) 


P = total precipitation during a time step, mm = P c + P p 

p C) P p = convective, large scale (uniform) precipitation rates during a time— step, 
respectively, mm 

The term PI^ can now be interpreted as an amount distribution which may be 
used to calculate throughfall rates and infiltration excess. 

Firstly, the direct throughfall term in SiB is defined as that fraction of the rainfall 
that falls through holes and gaps in the vegetation canopy. (Such holes and gaps are 
assumed to be uniformly spread over the grid area.) The definition of this direct fraction is 
given in Sellers et al. (1986), Sec. 7(b); here it is referred to as G; 0 < G < 1. 

Then, 

D d = GP (C.4) 

Dd = direct throughfall, mm 


PI/ \ — (P c a.c + P p a, p ) e* bx + (Pc Cc + Pp c p ) 



The amount of rainfall actually intercepted (though not necessarily retained) by the 
canopy is then given by 


= (l-G)PI, 


(C.5) 


The proportion of the grid area for which the canopy has intercepted enough 
rainfall to equal or exceed its saturation limit, x s , is then given by the solution to 


Y (xs) = S c -M c (C.6) 

S c = canopy storage limit, mm 

M c = intercepted water stored on canopy at the beginning of the time step, 
mm 
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Combining (C.6) with (C.3) we have 


x s 



log ■ 


S c — M c 

a(l-G) P 



(C.7) 


where a ~ (P c a c + Ppa p )/P and so forth. Reference to Fig. (C.2) will show the physical 
significance of x s ana how it is used to calculate the canopy drainage term, D c 


D c 


/ XS Y(x)dx-(S c -M c )x s 


= P(l-G)|(l-e- bx s) + cx s -(S c -M c )x s 


(C.8) 



S c 

M c 


Figure C.2: Dynamics of precipitation interception by a vegetation canopy in SiB. The 

amount of water already stored in the canopy prior to rainfall interception, 
Me, is considered to be uniformly distributed over the grid area (see hatched 
area). The integral of the water amount above M c represents to total 
amount of water intercepted by the canopy. x s is the proportion of the grid 
area where the intercepted rainfall plus the preexisting canopy water storage, 
M c , exceeds the canopy storage limit S c - All water above the 3 C limit drains 
off the canopy, all below it is added to the canopy interception store. 


Total throughfall to the ground, D, is then given by 
D — D<j + Dc 


(C.9) 
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The amount of water remaining on the canopy, i.e., available for direct evaporation 
is given by: 

M c = M c +P-D (C.10) 

t + 1 t 

For the next time step, this water is assumed to be uniformly spread over the 
vegetation canopy in the grid area which is clearly inconsistent with the representation 
shown above and the illustration shown in Fig. (C.2). However, one can assume that the 
resulting errors are small as the non— uniformity of intercepted precipitation will only 
determine the timing of evaporation loss unless convective storms are immensely 
concentrated in time and space. 

A similar calculation to that shown in (C.l) through (C.10) is performed for the 
ground cover so that an effective precipitation rate reaches the soil surface. With very 
heavy rainfall, the local infiltration capacity of the soil can be ‘beaten’ by the rainfall rate 
and (infiltration excess) overland flow can be generated. For the sake of simplicity, we 
assume the same area— amount relationship for the rainfall (in fact this would have been 
distorted by interception in the canopy, etc.). A soil hydraulic conductivity for the surface 
layer is calculated from Eq. 50 of Sellers et al. (1986) and this is used to calculate an 
equivalent depth of water that the soil can absorb during the time step. 

Pi = KAt (C.ll) 

Pi = amount of water that soil can absorb during the time step, mm 
K = soil hydraulic conductivity, mm s* 1 
At = time step duration, s 

The fraction of the grid area giving rise to infiltration excess is calculated by 
E S I( X1 ) = Pi (C.12) 

D s = amount of throughfall and vegetation drainage reaching the soil surface, 
mm 

The amount of infiltration excess is then simply given by 

Ri = /* Xl D s 1/ % dx — Pi (C.13) 

o*^ ' ' 

where the term inside the integral is similar to that in (C.8). 

When snow falls, Eqs. C.l to C.10 are used, but not C.ll to C. 13; that is, the 
snow can build up and fall off the vegetation, but accumulates on (and does not run off) 
the ground. 

This scheme is in place in all current SiB implementations in GCMs. 


62 



APPENDIX D 

INITIALIZATION OF GLOBAL SOIL MOISTURE FIELDS FOR THE 
SIMPLE BIOSPHERE MODEL (SiB) 


It was recognized in Sellers et al. (1986) that SiB would be sensitive to the soil 
moisture initialization as this can have a marked effect on the calculated 
evapotranspiration rate and hence on all the components of the surface energy balance. 
Unfortunately, the soil moisture climatologies currently used by GGMs are not directly 
compatible with SiB as the formulations used to derive them are drastically different from 
the biophysically— based calculations used in SiB. This appendix reviews a methodology 
that can be used to transform the existing soil moisture fields into SiB— compatible soil 
moisture fields; the method requires no other data than the existing soil moisture 
climatology and the distributed SiB parameter set. For a number of reasons the 
methodology has severe limitations. 


1. EXISTING SOIL MOISTURE CLIMATOLOGIES FOR GCM’S 

The fields of soil moisture content used for the initialization of GCMs originate 
from two principal sources. 


Mintz and Serafini (in prep.). Mintz and Serafini used a simple water 
budget equation set to calculate time varying global soil moisture fields on 
a 4° by 5‘ grid. 

The equations used were: 


a| = p-e , s < s„„ 

(D.l) 

E = /> ms E„ 

(D.2) 

/? us = l-exp -6 ' 8 ( S / S "“) 

(D.3) 


S = available soil moisture, mm 
Smax = maximum value of S; 150 mm 
P = rainfall rate, mm/day 
E = evapotranspiration rate, mm/day 
/?ms = Mintz and Serafini Beta-function which 
s r egulates the evapot ranspi ration r ate 
E p = potential evapotranspiration rate, mm/day. 
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The expression for /? Mg was determined by Mintz and Serafini by curve— fit to the 

data of Nappo (1975) which were based on observations of soil moisture depletion by 
individual maize plants growing in pots. The potential evapotranspiration rate, E p , was 
determined from the method of Thornthwaite (1948) and is principally a function of air 
temperature, T a . 

To calculate the time— varying fields of soil moisture, monthly mean values of P 
and T a were obtained from Jaeger (1976) and the NCAR data archive respectively and 
were used to integrate Eq. (D.l) lor each 4' by 5° grid area using a time step of one day. 
This procedure was repeated for 7 years of integration after which the difference between 
the values of S taken on the same calendar day of years 6 and 7 were less than 0.01 mm 
for every grid area. (The fields were initialized to zero on the first time step). These fields 
of soil moisture are currently used by the UCLA/GLA model and the GLA 4th-order 
model. In most GLA and UCLA/GLA GCM applications, the soil moisture fields are not 
updated prognostically. but from daily interpolations of the time— varying fields calculated 
by Mintz and Serafini (in preparation). 

(ii) Geophysical Fluid Dynamics Laboratory (GFDL) 

A similar scheme to that of Mintz and Serafini (in prep.) has been in 
use at GFDL and at the National Meteorological Center (NMC) for some 
years; see Gordon, (1986). 

A long term integration of the GFDL Zodiac GCM was used to 
generate fields of S via Eq. (D.l). In this case, the fields of P and E used 
in (D.l) are produced from the model’s own calculations where E is 
calculated by (D.2) with /? G given by 

h = inris < D - 4 ) 

where 

3 _ = GFDL ‘Beta— function’ 

u 

S m ax — 150 mm 


The potential evaporation rate, E p , is calculated as follows 
Ep = 


(D.5) 



la 


Pi Cp 


7 

A 

r 


potential evaporation rate, mm/day 
surface temperature, K 

saturated vapor pressure of the lowest model layer, mb 
aerodynamic resistance between the surface and the 
lowest model layer, s/m 

the density, specific heat of air, kg/m and J/kg/K, 
respectively 

psychrometric constant, mb/K 
latent heat of evaporation, J/kg 
number of seconds in a day. 
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The fields generated from the long-term GCM integration have been used as a 
standard soil moisture initialization for the GFDL model and the NMC medium range 
forecast (MRF) model. 

Recently, Willmott et al. (1985) performed another off-line analysis of the global 
soil moisture fields in a similar manner to that of Mintz and Serafini (in prep.). These 
fields (which are based on the /? Mg function) are currently used in the University of 

Maryland’s Center for Ocean— Land— Atmosphere version of the NMC MRF 86 model 
(Ctl-GCM) for initialization. For our purposes, however, we treat these fields as if they 
had been derived by the NMC bucket hydrological model (i.e., the 0 G function) as 

regardless of how the initial fields are derived, they impact the circulation through the 
hydrological model in place in the GCM. 

2. INCOMPATIBILITY OF EXISTING AND SiB MOISTURE CLIMATOLOGIES. 

The existing soil moisture initialization fields cannot be used directly by SiB for the 
following reasons: 

(i) The value of S max , the total amount of available soil moisture, is fixed at 150 
mm in the Mintz and Serafini and GFDL cases, but is variable in SiB depending 
on the biome. In SiB, for practical purposes the available soil moisture storage 
ranges from 60 mm in deserts to 415 mm in forested regions. 

(ii) The fractional wetness, W = S/S ma:c , as given by the existing fields, cannot be 
used for initializing SiB because the extremely low values of W predicted by 
(D.l) and (D.2) are physically unrealistic. Both of the Beta— functions, /? Bg 

and P Q , produce excessively high values of E/E p when W < 0.6, see Sellers 

(1987). The main reason for this is that vegetation cannot extract water from 
the soil when W drops below about 0.3. Most observations and simulations; 
for example, see Sellers and Dorman (1987), indicate that the evapotrans— 
piration rate of a vegetated surface is largely insensitive to W when W is 
greater than about 0.6, declines curvilinearly between 0.6 and about 0.3 and is 
near— zero thereafter. 

(iii) The values of E p as given by (D.4) and (D.5) will generally overestimate the 
stress-free evapotranspiration rate, see Sellers (1987) and Sellers et al. (1986). 
This effect is not likely to be so important as (i) and (ii). 

In conclusion it would be unwise to use the Mintz— Serafini or the GFDL— generated 
fields of soil moisture for SiB initialization directly. This applies to both the absolute 
values, S, and the wetness fraction, W. 

3. TRANSFORMING THE EXISTING SOIL MOISTURE FIELDS INTO 
SiB-COMPATIBLE SOIL MOISTURE FIELDS 

The transformation scheme used to generate SiB soil moisture fields is simple, 
robust and does not require any more data than the existing fields themselves and the SiB 
parameter set. The scheme is based on two representations of the soil moisture depletion 
process, one of which is assumed to be effective for near saturated conditions (Sec. 4.1) and 
the other for very dry conditions (Sec. 4.2). A simple interpolation formula is used to 
combine the two methods to determine intermediate soil moisture contents (Sec. 4.3). 
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3.1 Soil moisture depletion under wet conditions 

A given grid area is assumed to be saturated at time t = 0 and is then exposed to 
an evaporative demand, E p . After a specified time T, the grid area would have soil 

moisture contents of S Mg or S G , if the /? Mg or functions were applied, and a soil 

moisture content S g if SiB were used. The values of S Mg , and S g could be then 

considered transforms of each other. 

The process of integrating (D.l) over time is relatively straightforward. With the 
full version of SiB, however, the problem is more complex ana so the SiB evapo— 
transpiration formulation had to be considerably simplified so that it could be applied to 
the existing fields. To this end, the following assumptions were made. 

(i) The initial state of the Beta— function soil column is saturation; that is, 

S = Sjnax- The equivalent initial state assumed for the SiB soil column is the 

moisture content retained against gravitational drainage. This implies a soil 

moisture potential, ip ) of order —1, from which the appropriate soil moisture 
content can be estimated from 

ip = V’sW -0 (D.6) 

where 

ip = soil moisture potential fml 
ip s = soil moisture potential (m) when W = 1 
b = empirical coefficient 

Both ip s and b are dependent on soil type in SiB. For a given SiB biome, 
(D.6) may be solved for ip = — 1 to give an estimate of W 0 , the starting value 
of the SiB soil wetness when t = 0. 

(ii) The value of Ep used in the Beta— function calculations is taken to be roughly 
equivalent to the unstressed evapotranspiration rate from the equivalent SiB 
grid area. This is a less defensible assumption as firstly, the evapotranspiration 
rate from a SiB grid area is dependent on the vegetation type, and secondly, the 
value of E p used in the Beta— function is liable to be an over-estimate, see 
Sellers (1987). 

The SiB calculation of the surface energy balance can now be drastically simplified 
to yield a formulation similar to a Beta— function. To do this, some further assumptions 
are made: 

(iii) The ground and canopy surface temperatures are equal; T c = T gs 

(iv) One vegetation layer is dominant. 

(v) The leaf water potential is approximately equal to the soil moisture potential in 
the root zone, ip. 

(vi) Soil moisture stress, exercised through ip, is the only significant stress operating 
on the vegetation’s transpiration rate (see, Sellers et al. 1986). 

(vii) The surface resistances are large compared to the aerodynamic resistances. 

(viii) The air is relatively dry. 
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(D.T) 

(D.8) 


Incorporation of these assumptions into the SiB formulation yields 
E S " E p 

fs = Vf(« + (l-V)h s £aiio 



r surf 0 


*surf 


SiB— derived evapotranspiration rate, mm/day 
SiB— based Beta function 
dominant vegetation layer cover fraction 
soil air humidity, kg/kg. 

Soil surface resistance at saturation, s/m. 

Soil surface resistance, s/m 
leaf water potential stress factor. 


The leaf water potential stress factor, f (ijj), is given by 

« = (D.9) 

ipufa = leaf water potentials at which the stomata start to close 
and are completely closed, respectively, m. 

0 < f(tp) < 1 

Where ib is calculated from the empirical formulation of Clapp and Homberger (1978), see 
Eq. (D.6). 

The process of transforming the Beta— function fields is now fairly straightforward. 
The link between the Beta— function and SiB soil moisture contents is specified by the 
time, T, that each grid area is exposed to the same evaporative demand, E p . To do this, 
the precipitation rate in (D.l) is ignored which yields 

3T = n0Ep (D. 10) 


and 


where 


I 


MS,G , S 




dS 


MS , G,S 


“ — Ep / 


L dt 


(D-ll) 


So = 




starting soil moisture content, mm 
Smax for Beta— function calculations, mm 
1/B 


for SiB calculations, mm 


S = Soil moisture after time T, mm, where subscript refers to ms,g or SiB 


The SiB transform of the Mintz— Serafini or GFDL soil moisture fields is simply 
achieved by assuming E p is the same in both cases so that 




(D.12) 


In (D.12), the fields of S 110 _ and S 0 are known and the functions (3 r and 

MbjVJ Mb (j 

/?g have been specified in (D.3), (D.4), and (D.7) respectively. (D.12) can therefore be 
solved quite easily for S , the SiB soil moisture value. The left hand side of (D.12) is 
solved analytically to yield: 


and, 



Smax 


ln 


„6.8S us /S 


5 max 


e 6 • 8 -l 


0.75S ro ax 


In 


0.75S 


max 



4* 0.25Smaxj 


ax 



Sq < 0.75S ma x 
Sq ^ 0.75S m ax 


(D. 13) 


(DM) 


The integration of the right hand side of (D.12) is performed numerically to create 
a look— up table for each SiB biome. In the WETNESS program, soil moisture values from 
existing data sets are used to determine the left hand side of (D.12) for each grid area, 
using (D.13) and (D.14), and these are then used to obtain the corresponding values of S 

b 

from the look up table. 


3.2 Soil moisture depletion under dry conditions. 

In dry areas, it is assumed that the soil moisture content has reached a quasi- 
equilibrium state which is characterized by the value of the Beta— function, /? or 

Mb 


h- 


Under these conditions the soil moisture content for SiB, S D , is determined from the 
matching SiB soil moisture availability value as given by (D.6). 

,1/B 


S » = 


5max 




where ip is determined from the solution of 


Pc = A 


MS, 6 


(D.15) 


(D.16) 
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3.3 Interpolation of SiB soil moisture values for all wetness conditions. 

It is assumed that the methods outlined in Secs. 4.1 and 4.2 are appropriate to wet 
and dry conditions respectively. A single value of S is calculated simply by taking the 

b III 

minimum of S g and S D . 

S S,B = Mi "( S S>V (° 17 ) 


4. SUMMARY 

The method discussed above has a number of obvious shortcomings. It does serve, 
however, to give SiB— compatible soil moisture fields that will approximately reflect the 
evapotranspiration climatology of the region. The major features of the results from the 
procedure are as follows: 

(i) In mid-latitude zones where SiB has deep soils, the method will yield relatively 
high values of wetness fraction for SiB. This is because the soil moisture deficit 
given by (D.12) is calculated in absolute terms. 

(ii) Where the existing climatologies give very low values, i.e., W < 0.1, SiB will 
yield values of the order of W = 0.25, which is equivalent to unavailable soil 
moisture. 

(iii) The maximum value of the initial SiB wetness fraction will be about 
W = 0.80. 

All of the simplifications used to calculate the SiB fields ensure that a fairly high 
degree of uncertainty exists in these values. When a long integration with SiB is executed, 
the resultant fields of soil moisture will be used for future initializations. 
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APPENDIX E 

INITIALIZATION OF SOIL TEMPERATURE 


We have adopted a simple soil temperature initialization scheme for Ctl— GCM: 
the mean annual air surface temperature is used as the soil temperature at 500 cm depth, 
T g3 , the energy balance equation was then solved to obtain the ground surface 
temperatures for midday and night-time radiation conditions. These two surface 
temperatures were averaged to give the initial ground surface temperature, T s . The deep 
soil temperatures T gt and T g 2 (at 10 cm and 50 cm respectively) were initialized as 
follows 

T g i = 0.99 T s + 0.01 T g3 (E.l) 

T g2 = 0.9 T s + 0.1 T g 3 

In SiB-GCM, T g3 was used as the initial estimate of deep soil temperature, Td- SiB was 
then integrated for five time steps with noon time and zero solar radiation conditions while 
the atmospheric temperature and specific humidity were held constant. These two ground 
surface temperatures were averaged to give the initial surface temperature T s . The deep 
soil temperature Td was initialized as follows 

T d = 0.9 T s + 0.1 T g3 (E.2) 

As discussed in Sec. 4, this soil temperature initialization was found to be 
satisfactory and heat conduction in the soil was generally an insignificant part of the 
time-averaged surface energy balance. 
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flux and less latent heat flux over vegetated land than did Ctl-GCM and this 
resulted in the development of a much deeper daytime planetary boundary and reduced 
precipitation rates over the continents in SiB-GCM. In the summer simulation, the 
200 mb jet stream and the wind speed at 850 mb were slightly weakened in the SiB-GCMj 
relative to the Ctl-GCM results and equivalent analyses made from observations. 
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